-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathiceplume_calc.F
1037 lines (915 loc) · 42.4 KB
/
iceplume_calc.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include "ICEPLUME_OPTIONS.h"
CBOP
C !ROUTINE: ICEPLUME_CALC
C !INTERFACE:
SUBROUTINE ICEPLUME_CALC(
I myTime, myIter,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE ICEPLUME_CALC
C | o Send ambient conditions to plume model
C | o Calculate source/sink terms to parameterise movement of
C | water and tracers between vertical layers by resulting plume
C | o Calculate melt rates and tendencies due to melting of
C | ice front in none-plume locations
C |
C *==========================================================*
C | atn 09.Oct.2021: added notes
C | Because the thermodynamic tendencies code was copied from
C | pkg/icefront , the sign conventions are as follows:
C | Outward (pos) = leaving the ocean, inward (neg) = entering ocean.
C | Freshwater (kg/m2/s): inward (negative) freshwater flux implies
C | glacier melting due to outward (positive) heat flux (W/m^2)
C | kiki 27. October 2021: started to implement new plume option 6&7
C | frictionless plume following Wells and Worster, 2008 to simulate
C | purely melt water dricen plumes (see Jackson et al., 2020, GRL)
C | nad truncated line plume by Jackson et al., 2017
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "ICEPLUME.h"
#ifdef ALLOW_PTRACERS
#include "PTRACERS_PARAMS.h"
#include "PTRACERS_START.h"
#include "PTRACERS_FIELDS.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
_RL myTime
INTEGER myIter
INTEGER myThid
CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_ICEPLUME
C !LOCAL VARIABLES:
C == Local variables ==
C I,J,K,bi,bj :: loop indices
C msgBuf :: Informational/error message buffer
C sProf, tProf, uProf, vProf :: salt, pot. temperature and
C uVel and vVel profiles
C ptrProf :: ambient ptracer profile
C ptrPlumeCum :: cumulative quantity of ptracer in plume
C ptrPlume :: quantity of ptracer in plume outflow
C eps5 :: for thermodynamics (see pkg icefront)
C maxDepth :: vertical extent of domain (m)
C plumeAreaInCell :: surface area of plume in contact with ice in that cell (m^2)
C negSum, posSum :: sum of negative and positive contributions to the plume volume
C posNegRatio :: ratio of the above
C wVelCell :: vertical velocity component at cell centres
C hVelCell_tangential :: horizontal velocity component at cell centres
Catn vVelCell, uVelCell :: velocities at cell centres
C meanVel :: ice tangental velocity
catnC rho_0 :: average density of seawater
Catn lhfusion :: latent heat of fusion for solid to liquid phase, Joule
Catn rho_shelfice :: 917. kg/m3
C secInday :: number of seconds in a day
catn: introducing a new flag below to avoid confusion between where pos(neg)
catn: mask is NorthSouth or EastWest oriented glacier front
C GlacierFront_is_NorthSouth :: 1 for the glacier oriented NorthSouth
C :: 0 for EastWest, default is -9999
INTEGER bi, bj
INTEGER J, K, I
INTEGER GlacierFront_is_NorthSouth
_RL eps5
_RL plumeAreaInCell
_RL negSum, posSum, posNegRatio
_RL wVelCell, hVelCell_tangential, meanVel
_RL sw_temp, sw_ptmp
catn: There are too many loose rhos being defined. Here we make rho_0 = rhoConst
catn: Second, make rho_shelfice a readable param to avoid hardcoding.
catn _RL rho_0
catn _RL rho_shelfice
_RL secInday
external SW_TEMP
external SW_PTMP
catn PARAMETER(rho_0=1027.0D0)
PARAMETER(secInday=86400.0D0)
catn PARAMETER(rho_shelfice=917.0D0)
catn PARAMETER(lhfusion=333.55E+3)
catn: duplicate g[T,S]_iceplume from apply_forcing to see if the same
_RL gT_iceplume(1:sNx,1:sNy,Nr)
_RL gS_iceplume(1:sNx,1:sNy,Nr)
_RL gT_icefront(1:sNx,1:sNy,Nr)
_RL gS_icefront(1:sNx,1:sNy,Nr)
_RL mass_icefront(1:sNx,1:sNy,Nr)
_RL mass_iceplume(1:sNx,1:sNy,Nr)
#ifdef ALLOW_PTRACERS
INTEGER iTracer
_RL ptrPlume (PTRACERS_num)
_RL ptrPlumeCum (PTRACERS_num)
_RL ptrProf (Nr,PTRACERS_num)
#endif
catn: brought out of do-loop, factor for heat tendency
catn: mass2rUnit=1/rhoConst (set in data), e.g., 1/1029. [m3/kg]
catn: HeatCapacity_Cp = 3994. [J/kg/degC]
cat: so eps5 = 1/1029./3994. = 2.433193035422911e-07 [degC m3 / J]
eps5 = mass2rUnit/HeatCapacity_Cp
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_ENTER('ICEPLUME_CALC',myThid)
#endif
C -----------------------------------------
C Enter into loops
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
catn: initialize tendencies for plume
do k=1,Nr
do j=1,sNy
do i=1,sNx
gT_iceplume(i,j,k)=0. _d 0
gS_iceplume(i,j,k)=0. _d 0
gT_icefront(i,j,k)=0. _d 0
gS_icefront(i,j,k)=0. _d 0
mass_iceplume(i,j,k)=0. _d 0
mass_icefront(i,j,k)=0. _d 0
enddo
enddo
enddo
DO J = 1-OLy,sNy+OLy
DO I = 1-OLx,sNx+Olx
#ifdef ALLOW_PTRACERS
C Clear local plume ptracer variables
DO iTracer = 1,PTRACERS_num
ptrPlume(iTracer) = 0.D0
ptrPlumeCum(iTracer) = 0.D0
ENDDO
#endif /* ALLOW_PTRACERS */
catn initialize
GlacierFront_is_NorthSouth = -9999
C Check to see if there is ice in that cell. If not, skip to end.
IF ( plumeMask(I,J,bi,bj) .NE. 0 ) THEN
C Read from the plume mask which type of plume should be used in this cell.
C 1 = ice but no plume (melting only)
C 2 = sheet plume (Jenkins)
C 3 = half-conical plume (Morton/Slater)
C 4 = both sheet plume and half-conical plume (NOT YET IMPLEMENTED)
C 5 = detaching conical plume (Goldberg)
C 6 = buoyancy driven sublayer, Wells and Worster 2008 (Kiki)
C 7 = truncated line plume (Jacksom et al., 2017, GRL) (Kiki)
C POSITIVE values indicate ice front is orientated north-south
C NEGATIVE values indicate ice front is orientated east-west
catn: must use the sign as above to determine delta_[x,y]
C If there is subglacial discharge but no plume type defined, there will be no
C plume.
catn: this definition must go before the Q_sg
catn: moved this now to above calc of Q_sg.
catn: also, we MUST AVOID hardcoding delta_x and delta_y to [x,y] of
catn: model grid because it depends on the geometry of the plume, which
catn: is including in the sign of the plumeMask. So, since we're
catn: looping through I,J starting at the entry of this call, best to
catn: move def of delta_[x,y] to above when we're checking plumeMask.
catn: by default, it has been assummed that dyG is across the
catn: glacier and dxG is down the fjord toward open ocean. So
catn: we first define using that assumption:
C Cell resolution
IF ( plumeMask(I,J,bi,bj) .GT. 0 ) THEN
GlacierFront_is_NorthSouth = 1
dLnormal = dxG(I,J,bi,bj)
dLtangential = dyG(I,J,bi,bj)
ELSEIF ( plumeMask(I,J,bi,bj) .LT. 0 ) THEN
GlacierFront_is_NorthSouth = 0
catn: but now we check, if the plume geometry is such that dxG is across
catn: the glacier and dyG along the fjord toward open ocean:
dLtangential = dyG(I,J,bi,bj)
dLnormal = dxG(I,J,bi,bj)
ENDIF
IF ( plumeMask(I,J,bi,bj) .EQ. -1
& .OR. plumeMask(I,J,bi,bj) .EQ. 1) THEN
useSheetPlume = .FALSE.
useConePlume = .FALSE.
useDetachPlume = .FALSE.
useBuoyPlume = .FALSE.
useTruncPlume = .FALSE.
ELSEIF ( plumeMask(I,J,bi,bj) .EQ. -2
& .OR. plumeMask(I,J,bi,bj) .EQ. 2) THEN
useSheetPlume = .TRUE.
useConePlume = .FALSE.
useDetachPlume = .FALSE.
useBuoyPlume = .FALSE.
useTruncPlume = .FALSE.
ELSEIF ( plumeMask(I,J,bi,bj) .EQ. -3
& .OR. plumeMask(I,J,bi,bj) .EQ. 3) THEN
useSheetPlume = .FALSE.
useConePlume = .TRUE.
useDetachPlume = .FALSE.
useBuoyPlume = .FALSE.
useTruncPlume = .FALSE.
ELSEIF ( plumeMask(I,J,bi,bj) .EQ. -4
& .OR. plumeMask(I,J,bi,bj) .EQ. 4) THEN
useSheetPlume = .TRUE.
useConePlume = .TRUE.
useDetachPlume = .FALSE.
useBuoyPlume = .FALSE.
useTruncPlume = .FALSE.
ELSEIF ( plumeMask(I,J,bi,bj) .EQ. -5
& .OR. plumeMask(I,J,bi,bj) .EQ. 5) THEN
#ifdef ICEPLUME_ALLOW_DETACHED_PLUME
useSheetPlume = .FALSE.
useConePlume = .FALSE.
useDetachPlume = .TRUE.
useBuoyPlume = .FALSE.
useTruncPlume = .FALSE.
#else
WRITE(msgBuf,'(2A)')
& 'cannot use detaching plume without ',
& '#define ICEPLUME_ALLOW_DETACHED_PLUME'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R ICEPLUME_CALC'
#endif
ELSEIF ( plumemask(I,J,bi,bj) .EQ. -6
& .OR. plumemask(I,J,bi,bj) .EQ. 6) THEN
useSheetPlume = .FALSE.
useConePlume = .FALSE.
useDetachPlume = .FALSE.
useBuoyPlume = .TRUE.
useTruncPlume = .FALSE.
Ckiki can we change eq 7 to greater than some nbr and give Lp
Ckiki as input via the plume mask, Lp = width of trauncated line plume
Ckiki: also, remove options not implemented or used (4,5)
ELSEIF ( plumemask(I,J,bi,bj) .EQ. -7
& .OR. plumemask(I,J,bi,bj) .EQ. 7) THEN
useSheetPlume = .FALSE.
useConePlume = .FALSE.
useDetachPlume = .FALSE.
useBuoyPlume = .FALSE.
useTruncPlume = .TRUE.
ELSE
WRITE(msgBuf,'(2A)')
& 'Plume mask value must be between -5 and 5'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R ICEPLUME_CALC'
ENDIF
C Read in the subglacial discharge for this place and time
C (only the lowermost cell in column is read)
w_sg = runoffVel(I,J,bi,bj)
r_sg = runoffRad(I,J,bi,bj)
catn: there is a potential error for delta_y here because it is yet
catn initialized, and yet plays a role in Q_sg calc. Then after that we
catn define delta_y, and then compute (further down) volumeFlux based on
catn the redefined delta_y. So it has now been changed to dLtangential
IF ( useSheetPlume ) THEN
Q_sg = w_sg*r_sg*dLtangential
ELSEIF ( useConePlume ) THEN
Q_sg = w_sg*((pi*r_sg**2.)/2.)
ELSEIF ( useDetachPlume ) THEN
Q_sg = w_sg*((pi*r_sg**2.)/2.)
ELSEIF ( useBuoyPlume ) THEN
Q_sg = 1.0E-10
ELSEIF ( useTruncPlume ) THEN
ckiki: needs to be implemented using Lp, for now only placeholder
Q_sg = w_sg*r_sg
ELSE
Q_sg = 0.D0
ENDIF
c write(*,'(A,I10,1x,2I3,1X,F5.1,1X,F5.1,1X,F4.1,1X,3E15.7)')
c & 'iceplume_calc1: myIter,I,J,d[x,y],plumeMask,{w,r,Q}_sg: ',
c & myIter,I,J,dLnormal,dLtangential,plumeMask(I,J,bi,bj),
c & w_sg,r_sg,Q_sg
C Create variables with temperature, salinity
C and velocity profiles for that column
DO K = 1,Nr
C Tracers
catn prProf(k) = ABS(rC(k))*rho_0*9.81*1.0E-6 ! Pressure (dbar)
prProf(k) = ABS(rC(k))*rhoConst*gravity*1.0E-6 ! Pressure (dbar)
sProf(K) = salt(I,J,K,bi,bj) ! Salinity
ptProf(K) = theta(I,J,K,bi,bj) ! Potential Temperature
tProf(k) =
& SW_TEMP(sProf(k),ptProf(k),prProf(k),0. _d 0) ! Temperature
#ifdef ALLOW_PTRACERS
DO iTracer = 1,PTRACERS_num
ptrProf(k,iTracer) = pTracer(I,J,K,bi,bj,iTracer)
ENDDO
#endif /* ALLOW_PTRACERS */
C Velocities
vProf(k) = ABS(vVel(I,J,K,bi,bj)) ! v velocity
uProf(K) = ABS(uVel(I,J,K,bi,bj)) ! u Velocity
delta_z(k) = drF(K)
c write(*,'(A,I10,1x,I3,1X,F4.1,F9.2,1X,F5.2,1X,2F10.6,1X,2F8.2)')
c & 'iceplume_calc2: myIter,K,delta_z,{pr,s,pt,t,u,v}Prof:',
c & myIter,K,delta_z(k),prProf(k),sProf(k),ptProf(k),tProf(k),
c & uProf(k),vProf(k)
ENDDO
C Vertical extent of domain
maxDepth = rF(Nr+1)
C Depth of seabed at glacer front
iceDepth = R_low(I,J,bi,bj)
IF ( iceDepth .EQ. 0 ) THEN
WRITE(msgBuf,'(2A)')
& 'Plume specified in cell I = ', I, ', J = ', J,
& ', but depth of this cell = 0'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R ICEPLUME_CALC III'
endif
C Find grid layer at depth of iceplume
plumedepth = -ABS(plumedepth)
IF ( plumedepth .LT. iceDepth ) THEN
plumedepth = iceDepth
# ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('plumeD lt iceD',myThid)
# endif
endif
icedepthK = 0
DO K=1,Nr+1
IF ( rF(K) .EQ. plumedepth ) iceDepthK = K
ENDDO
C If a matching grid layer is not found, this may be because the bottom layer is a partial cell
C In this case, start in cell above partial cell
IF ( iceDepthK .EQ. 0 ) THEN
DO K=1,Nr+1
IF ( rF(K) .GT. plumedepth ) THEN
IF ( rF(K+1) .LT. plumedepth ) THEN
iceDepthK = K
ENDIF
ENDIF
ENDDO
ENDIF
# ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('icedepthkram4',myThid)
# endif
C If we still cannot find the bottom cell
IF ( iceDepthK .EQ. 0 ) THEN
# ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('icedepthk eq 0 WAAH',myThid)
# endif
WRITE(msgBuf,'(2A)')
& 'Unable to identify index of cell',
& 'at grounding line.',
& 'This may be because this is a partial cell.'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R ICEPLUME_CALC IV'
endif
C --- If there is subglacial outflow in that column, then parameterise plume ---
# ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('BEFORE Wsq gt 0 check',myThid)
# endif
IF ( Q_sg.GT.0 ) THEN
# ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_CALL('AFTER Wsq gt 0 check',myThid)
# endif
C This routine calculates T, S and W and r profiles for plume
CALL ICEPLUME_PLUME_MODEL (mythid)
C Calculate vertical plume volume flux...
DO k=1,Nr
C ... after checking to see if we are above the base of the ice face...
IF ( K .LT. iceDepthK ) THEN
C ... assuming specified plume horizontal extent (for sheet flow)...
IF ( useSheetPlume ) THEN
volFlux(k) = wProfPlume(k)*rProfPlume(k)*dLtangential
C ... or assuming half-conical form
ELSEIF ( useConePlume ) THEN
volFlux(k)=pi*(rProfPlume(k)**twoRL)*wProfPlume(k)/twoRL
C ... assuming detached-conical form
ELSEIF ( useDetachPlume ) THEN
volFlux(k)=rProfPlume(k)*wProfPlume(k)
C ... assuming buoyancy driven sublayer plume (same as sheet plume)
ELSEIF ( useBuoyPlume ) THEN
volFlux(k) = wProfPlume(k)*rProfPlume(k)*dLtangential
C ... assuming truncated line plume
ckiki: needs to be implemented using Lp, for now placeholder
ELSEIF ( useTruncPlume ) THEN
volFlux(k) = wProfPlume(k)*rProfPlume(k)*dLtangential
ENDIF
ELSE
volFlux(k) = 0.D0
ENDIF
c write(*,'(A,I10,1x,I5,1x,F14.6)')
c & 'iceplume_calc4: myIter,k,volFlux: ',myIter,k,volFlux(k)
ENDDO
C A couple of corrections:
C - even if plume is still buoyant, it cannot flow through the fjord surface
volFlux(1) = 0.D0
C - the initial volume flux is equal to runoff [m3/s]
volflux(iceDepthK) = Q_sg
c write(*,'(A,I10,1x,3I6,1X,2E15.7)')
c & 'iceplume_calc5: myIter,I,J,iceDepthK,volFlux(1,iceDepthK):'
c & ,myIter,I,J,iceDepthK,volFlux(1),volFlux(iceDepthK)
C Calculate volume flux differential to give entrainment / extrainment, [m3/s]
C First clear volfluxdiff
DO K = 1,Nr
volfluxdiff(K) = 0.D0
ENDDO
DO k=1,iceDepthK-1
volFluxDiff(k) = volFlux(k+1) - volFlux(k)
c write(*,'(A,I10,1x,I6,1x,E15.7)')
c & 'iceplume_calc6: myIter,k,volFluxDiff: ',
c & myIter,k,volFluxDiff(k)
ENDDO
catn: The loop below indicates that when conserveMass is true, then we scale
catn: the volume flux (in the plume) total to reduce the effect of Q_sg into
catn: the domain. An external test shows that before entering this loop, we
catn: have a non-zero sum(volFluxDiff), but after entering and getting
catn: scaled, sum(volFluxDiff) = 0 (within precision), and any posSum at
catn: any particular klev is scaled while all negSum were untouched.
IF ( conserveMass ) THEN
C Scale output to compensate for entrainment lost in expanding of output layer
C i.e. so that there is no net flow over boundary
catn: sign convention: if volFluxDiff(k) is less than 1, there is a net
catn: divergence == volume loss at the klev, and will contrib to negSum
negSum = 0.D0
posSum = 0.D0
DO K = 1,Nr
IF ( volFluxDiff(K) .LT. 0 ) THEN
negSum = volFluxDiff(K) + negSum
ELSE
posSum = volFluxDiff(K) + posSum
ENDIF
ENDDO
catn: the check of posSum not eq 0 implies there is at least ONE klev where
catn: there is a net CONVERGENCE into the cell in one or more klev. So we
catn: take the ratio of -negSum/posSum to check for the total columm
catn: what the net convergence is. If this ratio is one, negSum=posSum
catn: and there is no net water excess. If posSum>negSum, there is
catn: excess water input that needs to be taken care of under the
catn: scenario where ocean volume cannot be changed. So take the
catn: example negSum = 4 and posSum = 5, ratio is 4/5. This ratio will
catn: now be mult across all volFluxDiff>0, such that after scaling the
catn: net column convergence is zero
IF ( posSum .NE. 0 ) THEN
posNegRatio = -negSum / posSum
DO K = 1,Nr
IF ( volFluxDiff(K) .GT. 0 )
& volFluxDiff(K) = volFluxDiff(K) * posNegRatio
ENDDO
ENDIF
ENDIF
#ifdef ALLOW_PTRACERS
C Add up total sum of each tracer in plume
DO K=1,iceDepthK-1
IF (volFLuxDiff(k) .LT. 0. ) THEN
DO iTracer = 1,PTRACERS_num
ptrPlumeCum(iTracer)
& = ptrPlumeCum(iTracer)
& +(-volFluxDiff(k)*ptrProf(k,iTracer))
ENDDO
ENDIF
ENDDO
C Add ptracers in runoff
IF ( useInputPtracers ) THEN
DO iTracer = 1,PTRACERS_num
IF (ptracerMask(I,J,iTracer,bi,bj) .NE. 0 ) THEN
ptrPlumeCum(iTracer) =
& ptrPlumeCum(iTracer) +
& ptracerMask(I,J,iTracer,bi,bj) * ! ptracerMask is now a nx by ny by n_ptracers matrix
& volFlux(iceDepthK)
ENDIF
ENDDO
ENDIF
C Calculate concentration of tracer in outflow
DO K=1,iceDepthK-1
IF (volFluxDiff(k) .GT. 0. ) THEN
DO iTracer = 1,PTRACERS_num
ptrPlume(iTracer)
& = ptrPlumeCum(iTracer) / volFluxDiff(k)
ENDDO
ENDIF
ENDDO
#endif /* ALLOW_PTRACERS */
ELSE ! ( Q_sg .EQ. 0 )
C If no subglacial output, then there is no plume
DO k = 1,Nr
rProfPlume(K) = 0.D0
wProfPlume(K) = 0.D0
tProfPlume(K) = 0.D0
sProfPlume(K) = 0.D0
aProfPlume(K) = 0.D0
mIntProfPlume(K) = 0.D0
#ifdef ICEPLUME_ALLOW_DETACHED_PLUME
thetaProfPlume(k) = 0.D0
distanceProfPlume(k) = 0.D0
#endif
ENDDO
ENDIF ! ( Q_sg.NE.0 ) THEN
C Send outputs to 3D grid for diagnostics
IF ( usePlumeDiagnostics ) THEN
DO K = 1,Nr
rProfPlume3D(I,J,K,bi,bj) = rProfPlume(k)
wProfPlume3D(I,J,K,bi,bj) = wProfPlume(k)
tProfPlume3D(I,J,K,bi,bj) = tProfPlume(k)
sProfPlume3D(I,J,K,bi,bj) = sProfPlume(k)
#ifdef ICEPLUME_ALLOW_DETACHED_PLUME
thetaProfPlume3D
& (I,J,K,bi,bj) = thetaProfPlume(k)
distanceProfPlume3D
& (I,J,K,bi,bj) = distanceProfPlume(k)
#endif
ENDDO
ENDIF
C-------- Calculate melt rates ----------------------------
DO K = 1,Nr
C Check to see if we are above the sea bed
IF ( K .GE. iceDepthK ) THEN
C If not then there is no melting
mProfAv(k) = 0.D0
mProfPlume(k) = 0.D0
mProf(k) = 0.D0
ELSE !k le iceDepthK
C If there is a plume in that cell, then need to calculate plume melt rate [m/time]
C distinct to background melt rate. Plume melt rate is already encorporated in
C the plrume model, and taken into account in the temperature and salinity of the
C plume outflow. It is useful though to have it available as a diagnostic.
plumeAreaInCell = 0.0
catn: note that mProfPlume below is not impacted by scaling of volFluxDiff
IF ( ( Q_sg .NE. 0 ) .AND.
& (useConePlume .OR. useSheetPlume.OR.useDetachPlume .OR.
& useBuoyPlume) )
& THEN
catn: aProfPlume(k): integrated contact area, calc in iceplume_plume_model [m2]
plumeAreaInCell = aProfPlume(k) - aProfPlume(k+1)
catn: if the area above is bigger than below, implying plume area expands
IF (plumeAreaInCell .gt. 0.0) then
catn: mIntProfPlume: integrated melt rate [m/day]
mProfPlume(k) =(mIntProfPlume(k)-mIntProfPlume(k+1))/
& plumeAreaInCell
& * secInday
ELSE
mProfPlume (k) = 0.0
ENDIF
c write(*,'(A,A,I10,1x,I4,1x,6E15.7)')
c & 'iceplume_calc7: myIter,k,plumeAreaInCell,mProfPlume(k),'
c & ,'mIntProfPlume(k+1,k),[t,s]ProfPlume(k): ',myIter,k,
c & plumeAreaInCell,mProfPlume(k),mIntProfPlume(k+1),
c & mIntProfPlume(k),tProfPlume(k),sProfPlume(k)
ELSE !Q_sg eq 0
C If there is no plume in that cell, set plume melt rate to zero
mProfPlume(k) = 0.D0
ENDIF
C Calculate the background melt rate (i.e. not generated by plumes). This will
C then be used to update the temperature and salinity in the adjacent cells.
C Velocities are calculated at cell faces - find averages for cell centres.
C Does not include velocity perpendicular to ice - this differs depending on
C orientation of ice front
catn: compute only one horz tangential vel based on ice front orientation:
catn: if North-South (East-West), we use vvel (uvel)
IF ( GlacierFront_is_NorthSouth .EQ. 1 ) THEN
hVelCell_tangential = (ABS(vVel(I,J,K,bi,bj))
& +ABS(vVel(I,J+1,K,bi,bj))) / twoRL
ELSEIF ( GlacierFront_is_NorthSouth .EQ. 0 ) THEN
hVelCell_tangential = (ABS(uVel(I,J,K,bi,bj))
& +ABS(uVel(I+1,J,K,bi,bj))) / twoRL
ELSE
catn: for adjoint
hVelCell_tangential = 0. _d 0
ENDIF
IF ( K .LT. Nr ) THEN
wVelCell = (ABS(wVel(I,J,K,bi,bj))
& +ABS(wVel(I,J,K+1,bi,bj))) / twoRL
ELSE
wVelCell = ABS(wVel(I,J,K,bi,bj)) / twoRL
ENDIF
catn: From above:
C POSITIVE values indicate ice front is orientated north-south
C NEGATIVE values indicate ice front is orientated east-west
catn: So there is a discrepancy between above and below. The importance
catn: is to settle on one convention and remain consistent. For now, we
catn: stay with convention established above, now improved with a flag.
catn IF ( plumeMask(I,J,bi,bj) .LT. 0 ) THEN
catn-C Negative mask values indicate east-west ice front orientation
catn IF ( GlacierFront_is_Northsouth .EQ. 0 ) THEN
catn meanVel = ((wVelCell**twoRL)+(uVelCell**twoRL))**HalfRL
catn- ELSEIF ( plumeMask(I,J,bi,bj) .GT. 0 ) THEN
catn-C Positive mask values indicate north-south ice front orientation
catn ELSEIF ( GlacierFront_is_NorthSouth .EQ. 1 ) THEN
catn meanVel = ((wVelCell**twoRL)+(vVelCell**twoRL))**HalfRL
catn ENDIF
catn: From Tom: what we want here is meanVel tangential to the ice front.
meanVel = ((wVelCell**twoRL)
& +(hVelCell_tangential**twoRL))**HalfRL
CALL ICEPLUME_MELTRATE(
I tProf(k),sProf(k),meanVel,rC(k),
O mProf(k) )
c write(*,'(A,I10,1x,3I5,1X,3E15.7,1x,F8.2,1x,E15.7)')
c & 'iceplume_calc8: myIter,I,J,K,{t,s}Prof,meanVel,rC,mProf:',
c & myIter,I,J,K,tProf(k),sProf(k),meanVel,rC(k),mProf(k)
C Get average melt rate. This is useful for visualling melt patterns and
C assessing overall melt rate of glacier. Unit: [m/time]
C the following should apply to both conical and sheet plume models
IF ( ( Q_sg .NE. 0 ) ) THEN
plumeAreaInCell = aProfPlume(k) - aProfPlume(k+1)
catn: this is quite troublesome here that delta_y kept entering as if
catn: this geometry is strictly applied for only ice front aligning
catn: along u-direction such that the cross-sectional dlength is delta_y
catn: Now we change delta_y to dLtangential
IF ( plumeAreaInCell .LE. dLtangential*delta_z(k) ) THEN
IF ( plumeAreaInCell .LE. 0 ) THEN
C If there is no plume in cell, then the melt rate is equal to the background melt rate.
mprofAv(k) = mProf(K)
ELSE
C If there is a plume in cell, calculate average melt rate
catn: again the trouble with delta_y below;
catn: We now change delta_y to dLtangential later
mProfAv(k) = (mProfPlume(k)*plumeAreaInCell
& +mProf(k)*(dLtangential*delta_z(k)-plumeAreaInCell))
& / (dLtangential * delta_z(k))
C Scale down background melt rate to account for area occupied by plume
C (necessary so that tendency terms aren't over estimated)
mProf(k) = mProf(k)*(1-plumeAreaInCell/
& (dLtangential*delta_z(k)))
ENDIF
ELSE
C if the plume contact area is larger than the cell area, we assume there is
C no background melting
mProfAv(k) = mProfPlume(k)*plumeAreaInCell /
& (dLtangential*delta_z(k))
mProf(k) = 0.
ENDIF
ELSE ! Q_sg eq 0, not coneplume or sheet plume
C If it is not a plume cell, then no plume melting.
mProfPlume(k) = 0.D0
mProfAv(k) = mProf(k)
ENDIF ! plume type
c write(*,'(A,A,I10,1x,3I5,1x,2F6.1,1x,4E15.7)')
c & 'iceplume_calc9: myIter,I,J,K,delta{y,z},plumeAreaInCell,',
c & 'mProfPlume,mProf,mProfAv: ',myIter,I,J,K,dLtangential,
c & delta_z(k),plumeAreaInCell,mProfPlume(k),mProf(k),mProfAv(k)
ENDIF ! are we above sea bed?
C Send outputs to 3D grid for diagnostics
catn: Since we're saving these 3D diags, should know how to calc offline
catn: their contributions to addMass3D(plume,bg)
IF ( usePlumeDiagnostics ) THEN
mProfPlume3D(I,J,K,bi,bj) = mProfPlume(k)
mProfAv3D(I,J,K,bi,bj) = mProfAv(k)
ENDIF
C ------------Tendencies------------------------------
C These are applied in cells where there is no plume. The idea is that in these areas there are likely to be
C local subgrid convection cells. As such, it is most realistic to apply changes in T and S to ice edge cell.
C Where there is a plume, products of melt are transported upwards so no local changes applied.
C The thermodynamics in this section are taken from pkg/icefront
C Tendencies are applied in S/R EXTERNAL_FORCING
C To convert from melt rate (m d^-1) to freshwater flux (kg m^-2 d^-1)
catn: m/day * kg/m3 = kg / day / m^2. where does 94.22 come from? I also think the
catn: quoted unit above is incorrect and fwflux is in kg/m2/s, same as in icefront.
catn: rhoshelfice/secInday = 917./(24*3600) = 0.0106 = 1/94.22
catn: sign convention from pkg/icefront: inward (negative, into the ocean) fwFlux
catn: implies glacier melting due to outward (positive, leaving the ocean) heatflux
catn: The last step we need to take care of is that FwFlux(k) is actually computed
catn: using rhoShelfIce, but outside of this addMass* is converted to Eta using
catn: rUnit2mass = rhoConst. How to reconcile this? Do we mult FwFlux back by
catn: rhoShelfIce to convert to m3/s, then rUnit2mass? If we look at how eccov4r5
catn: pkg/shelfice is treating fwflux from background melt, they use rUnit2mass and
catn: NOT rhoShelfIce. In light of this, we REMOVE rhoShelfIce from the calculation
catn: of FWflux and use rUnit2Mass=rhoConst.
catn FwFlux(k) = -mProf(k)/94.22
catn FwFlux(k) = -mProf(k)*rhoShelfIce/secInday
FwFlux(k) = -mProf(k)*rUnit2Mass/secInday
catn: note above: if we take it from merged pkg/shelfice and pkg/icefront, the
catn: melt is computed following Holland and Jenkin 1991 as something like:
catn: mProf=1/SHELFICElatentHeat*(
catn: SHELFICEheatCapacity_Cp*SHELFICEkappa/iceConductionDistance*(Tfreeze-thetaIceConduction)
catn: - rhoConst/rhoShelfice*HeatCapacity_Cp*gammaT*(Tinsitu-Tfreeze) ) [m/s]
catn: This means if we looking into iceplume_meltrate.F, we should see something like:
catn: SHELFICEheatCapacity_Cp (2000 J/kg/K), SHELFICEkappa (molecular thermal conductivity,
catn: 1.54e-6 m2/s), iceConductionDistance (~100m), thetaIceConduction (-20. degC) ?
catn: Lastly, to convert this FwFlux (kg/m2/s) to addMass3Dbg [kg/s], mult by the vertical
catn: area (dLtangential*dz). Note we have already taken care at this point to ensure
catn: dLtangential is along the glacier front. In addition, we've already taken care to
catn: reduce this addmass3dbg to zero if the whole cell is plume.
catn addMass3Dbg(I,j,k,bi,bj)=Fwflux3D(I,J,k,bi,bj)*dLtangential*delta_z(k)*rUnit2mass
catn: pay attention to sign! addMass is minus FwFlux!! (following pkg/shelfice of eccov4r5)
addMass3Dbg(I,J,k,bi,bj)=-FwFlux(k)*dLtangential*delta_z(k)
C Heat required to melt that much ice (W m^-2)
catn: 333.55 x 10^3 is the Latent heat of fusion: When melting 1kg of ice, 333.55 kJ of
catn: energy is absorbed with no temperature change. While melting, the heat energy
catn: needed to change a substance from solid to liquid at atm pressure is the latent
catn: heat of fusion. The liquid phase has a higher internal energy than the solid phase.
catn: From statement above, in an ocean-ice system, when there is an dheff ice melt, implying
catn: -dheff*rhoi/rhosw=fwflux (negative) to the ocean, the ocean first needs to put energy
catn: into raising the ice temperature to freezing, then put more energy into solid->liquid
catn: phase change. Thus, from the ocean perspective, it loses dheat to raise ice temp plus
catn: dheat_Latent convert solid->liquid. So a +dheff gives -fwflx (neg, inward, into ocean),
catn: the ice system gains energy and the ocean loses energy and cools down
catn: (outward, postive=leaving ocean). dheat_lost_by_ocean is opposite sign of fwflux.
catn heatflux(k) = -FwFlux(k)*333.55E+3
heatflux(k) = -FwFlux(k)*ICEFRONTlatentHeat
C Create local (no overlap) arrays of heat and freshwater flux from background melting
catn: more problematic assumption about I vs J and geometry of glacier front
IF (GlacierFront_is_NorthSouth .eq. 1) THEN
IF ( ( J .GT. 0 ) .AND. ( J .LT. sNy+1 ) ) THEN
Fwflux3D(I,J,k,bi,bj) = FwFlux(k)
HeatFlux3D(I,J,k,bi,bj) = HeatFlux(k)
ENDIF
ELSEIF (GlacierFront_is_NorthSouth .eq. 0) THEN
IF ( ( I .GT. 0 ) .AND. ( I .LT. sNx+1 ) ) THEN
Fwflux3D(I,J,k,bi,bj) = FwFlux(k)
HeatFlux3D(I,J,k,bi,bj) = HeatFlux(k)
ENDIF
ENDIF
C Compute tendencies (as for pkg/icefront)
catn: In ini_parms.F: mass2rUnit = recip_rhoConst, e.g., 1/999.8=1.0002e-3, HeatCapacity_Cp=3994
catn: Sign convention: using pkg/icefront as a guide: ocean will take the negative of the
catn: heatflux and positive of fwflux above, i.e., opposite sign of what was computed above.
icefront_TendT(I,J,K,bi,bj) =
& - HeatFlux3D(I,J,K,bi,bj)*eps5
icefront_TendS(I,J,K,bi,bj) =
& FWFlux3D(I,J,K,bi,bj)*
& mass2rUnit * sProf(k)
C Scale by icefrontlength, which is the ratio of the horizontal length
C of the ice front in each model grid cell divided by the grid cell area.
C (icefrontlength = dy / dxdy = 1 / dx)
catn: more problem with assuming geometry delta_x and delta_y
catn: solved now by rename to dLtangential and dLnormal
catn: scale factor: dLtangential/(dLnormal*dLtantengial) = 1/dLnormal
icefront_TendT(I,j,K,bi,bj) =
& icefront_TendT(I,j,K,bi,bj)
& * 1./dLnormal
icefront_TendS(I,j,K,bi,bj) =
& icefront_TendS(I,j,K,bi,bj)
& * 1./dLnormal
c write(*,'(A,A,I10,1x,I6,1x,5E15.7)')
c & 'iceplume_calc10: myIter,K,FwFlux,heatflux,',
c & 'icefront_Tend[T,S],addMass3Dbg: ',myIter,K,FwFlux(k),
c & heatflux(k),icefront_TendT(I,J,K,bi,bj),
c & icefront_TendS(I,J,K,bi,bj),addMass3Dbg(I,J,K,bi,bj)
ENDDO !the big k-loop
catn: Note: currently we're not taking care of the mass added from FwFlux (background melt).
catn: To fully take into account all fw added to the ocean, we will need to either add it
catn: to addMass array, or create yet another 3D field specific to iceplume where we store
catn: it there. The latter, though not desirable as we're carring more 3d fields around,
catn: is very good for diagnostic purpose, as we can for sure compute the ratio of
catn: background melt vs plume-induced entrainment vs Qsg
C The plume transport is undertaken using the addMass terms.
C addMass terms for volume (kg/s), pot. temperature, salinity and ptracers are input
C into the correct locations in 3D arrays
catn: need to first compute addMass specific to plume and store it independently IFF we
catn: want to use separate processes to treat Tendencies at different places in the code
IF ( Q_sg.NE.0. ) THEN
C Find temperature and salinity of plume outflow
DO k = 1,Nr !second kloop
catn: volFluxDiff entering here: when there's net convergence at this k
IF ( volFluxDiff(k) .GT. 0. ) THEN
temp_AddMass3D(I,J,K,bi,bj) = ! convert to potential temp
& SW_PTMP(sProfPlume(k),tProfPlume(k),prProf(k),0. _d 0)
salt_AddMass3D(I,J,K,bi,bj) = sProfPlume(k)
#ifdef ALLOW_PTRACERS
DO iTracer = 1,PTRACERS_num
ptr_AddMass3D(I,J,K,bi,bj,iTracer)
& = ptrPlume(iTracer)
ENDDO
#endif /* ALLOW_PTRACERS */
catn
ELSE
temp_AddMass3D(I,J,K,bi,bj) = ptProf(k)
salt_AddMass3D(I,J,K,bi,bj) = sProf(k)
#ifdef ALLOW_PTRACERS
DO iTracer = 1,PTRACERS_num
ptr_AddMass3D(I,J,K,bi,bj,iTracer)
& = ptrProf(k,iTracer)
ENDDO
#endif /* ALLOW_PTRACERS */
ENDIF
catn: now also filling out array for background melt, although not needed for
catn: tendency because we already take care of it in here with tend[T,S]
C Convert m3/s into kg/s
catn: volFluxDiff is now converted to addMass, so already taken care of massflux into ocean at depth
catn: Need to avoid hardcoded 1000 (missing decimal), will use rhoConstFresh
catn: Note: in integr_continuity, we convert addMass to volume using rhoConst=1029, so we'll need to
catn: use that here to avoid a difference in "mass" effect between what's computed in here and outside
catn addMass(I,j,k,bi,bj) = volFLuxDiff(k)*1000
catn addMass(I,j,k,bi,bj) = volFLuxDiff(k)*rhoConstFresh
catn: rUnit2mass = rhoConst, this is how mitgcm converts between volume and mass
addMass3Dplume(I,j,k,bi,bj) = volFLuxDiff(k)*rUnit2mass
c write(*,'(A,A,I10,1x,3I5,1x,F12.7,1x,F12.7,1x,E15.7)')
c & 'iceplume_calc11: myIter,I,J,K,[temp,salt]_AddMass3D,',
c & 'addMass3Dplume: ',myIter,I,J,K,temp_AddMass3D(I,J,K,bi,bj)
c & ,salt_AddMass3D(I,J,K,bi,bj),addMass3Dplume(I,j,k,bi,bj)
ENDDO !kloop
ELSE !Q_sg eq 0
DO K = 1,Nr
temp_AddMass3D(I,J,K,bi,bj) = ptProf(k)
salt_AddMass3D(I,J,K,bi,bj) = sProf(k)
catn: cannot reset a global field in here, need to make process-specfic field
catn addMass(I,j,k,bi,bj) = 0.D0
addMass3Dplume(I,j,k,bi,bj) = 0.D0
c write(*,'(A,A,I10,1x,3I5,1x,F12.7,1x,F12.7,1x,E15.7)')
c & 'iceplume_calc11: myIter,I,J,K,[temp,salt]_AddMass3D,',
c & 'addMass3Dplume: ',myIter,I,J,K,temp_AddMass3D(I,J,K,bi,bj)
c & ,salt_AddMass3D(I,J,K,bi,bj),addMass3Dplume(I,j,k,bi,bj)
ENDDO
ENDIF
ENDIF ! plumeMask .NE. 0
C DO J loop
ENDDO
C DO I loop
ENDDO
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
C Save plume values
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics .AND. usePlumeDiagnostics ) THEN
C Transfer to local (no bi,bj indices) and interior only arrays
C (can't seem to get the diagnostics to output properly with the
C other options!)
DO K = 1,Nr
DO J = 1,sNy
DO I = 1,sNx
wProfPlume3dLocal(I,J,K) =
& wProfPlume3D(I,J,K,bi,bj)
tProfPlume3dLocal(I,J,K) =
& tProfPlume3D(I,J,K,bi,bj)
sProfPlume3dLocal(I,J,K) =
& sProfPlume3D(I,J,K,bi,bj)
rProfPlume3dLocal(I,J,K) =
& rProfPlume3D(I,J,K,bi,bj)
mProfPlume3dLocal(I,J,K) =
& mProfPlume3D(I,J,K,bi,bj)
mProfAv3dLocal(I,J,K) =
& mProfAv3D(I,J,K,bi,bj)
#ifdef ICEPLUME_ALLOW_DETACHED_PLUME
thetaProfPlume3DLocal(I,J,K) =
& thetaProfPlume3D(I,J,K,bi,bj)
distanceProfPlume3dLocal(I,J,K) =
& distanceProfPlume3D(I,J,K,bi,bj)
#endif
mass_icefront(i,j,k)=addMass3Dbg(i,j,k,bi,bj)
mass_iceplume(i,j,k)=addMass3Dplume(i,j,k,bi,bj)
gT_icefront(i,j,k)=icefront_tendT(i,j,k,bi,bj)
gS_icefront(i,j,k)=icefront_tendS(i,j,k,bi,bj)
IF ( selectAddFluid.NE.0) THEN
IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
& .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
gT_iceplume(i,j,k)=addMass3Dplume(i,j,k,bi,bj)*mass2rUnit
& *( temp_addMass3D(I,J,k,bi,bj) - theta(i,j,k,bi,bj) )
gT_iceplume(i,j,k)=addMass3Dplume(i,j,k,bi,bj)*mass2rUnit
gT_iceplume(i,j,k)=addMass3Dplume(i,j,k,bi,bj)*mass2rUnit
& *( temp_addMass3D(I,J,k,bi,bj) - tRef(k) )
& *recip_rA(i,j,bi,bj)
& *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
gS_iceplume(i,j,k)=addMass3Dplume(i,j,k,bi,bj)*mass2rUnit
& *( salt_addMass3D(I,J,k,bi,bj) - SRef(k) )
& *recip_rA(i,j,bi,bj)
& *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
ENDIF
ENDIF
ENDDO !I
ENDDO !J
ENDDO !K
C Output diagnostics
DO K=1,Nr
C Here, the first 'k' is the layer in the output field in which to save the data
C and the second k is layer from which this data is taken in the original field
CALL DIAGNOSTICS_FILL_RS(wProfPlume3DLocal,'icefrntW',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(tProfPlume3DLocal,'icefrntT',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(sProfPlume3DLocal,'icefrntS',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(rProfPlume3DLocal,'icefrntR',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(mProfPlume3DLocal,'icefrntM',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(mProfAv3DLocal,'icefrntA',
& k,k,3,bi,bj,myThid)
#ifdef ICEPLUME_ALLOW_DETACHED_PLUME
CALL DIAGNOSTICS_FILL_RS(thetaProfPlume3DLocal,'PlumAngl',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(distanceProfPlume3dLocal,'PlumDist',
& k,k,3,bi,bj,myThid)
#endif
CALL DIAGNOSTICS_FILL_RS(gT_icefront,'IP_gTbg ',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(gS_icefront,'IP_gSbg ',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(gT_iceplume,'IP_gTplm',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(gS_iceplume,'IP_gSplm',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(mass_iceplume,'IPmasspl',
& k,k,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(mass_icefront,'IPmassbg',
& k,k,3,bi,bj,myThid)
ENDDO
C Clear local arrays otherwise results replicate on other tiles
DO I = 1,sNx
DO J = 1,sNy