New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
traicw.F90 in branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/NERC/dev_r5589_marine_glacier_termini/NEMOGCM/NEMO/OPA_SRC/TRA/traicw.F90 @ 5605

Last change on this file since 5605 was 5605, checked in by mathiot, 9 years ago

Marine glacier termini: initial commit
File size: 33.4 KB
Line 
1MODULE traicw
2   !!==============================================================================
3   !!                       ***  MODULE  trasbc  ***
4   !! Ocean active tracers:  surface boundary condition
5   !!==============================================================================
6   !! History : 02/2015  P. Mathiot : original code
7   !!----------------------------------------------------------------------
8   !!----------------------------------------------------------------------
9   !!   tra_icw      : parametrisation of ice wall processes
10   !!----------------------------------------------------------------------
11   USE oce             ! ocean dynamics and active tracers
12   USE dom_oce         ! ocean space domain variables
13   USE in_out_manager  ! I/O manager
14   USE fldread         ! read input field at current time step
15   USE iom
16   USE wrk_nemo        ! Memory Allocation
17   USE timing          ! Timing
18   USE eosbn2
19
20   IMPLICIT NONE
21   PRIVATE
22
23   PUBLIC   tra_icw    ! routine called by step.F90
24   PUBLIC   tra_icw_alloc ! routine call in sbcmod module
25   PUBLIC   tra_icw_init  ! (PUBLIC for TAM)
26   PUBLIC   div_icw   ! routine called in sshwzv module
27
28   LOGICAL, PUBLIC  :: ln_plumes
29   LOGICAL, PUBLIC  :: ln_ovt
30   LOGICAL, PUBLIC  :: ln_traicw
31
32   CHARACTER(len=100), PUBLIC ::   cn_dir          !: Root directory for location of icw files
33
34   TYPE(FLD_N)       , PUBLIC ::   sn_qmelt
35   TYPE(FLD_N)       , PUBLIC ::   sn_qsubg
36   TYPE(FLD_N)       , PUBLIC ::   sn_ovt
37   TYPE(FLD_N)       , PUBLIC ::   sn_ztop
38   TYPE(FLD_N)       , PUBLIC ::   sn_zovt
39   TYPE(FLD_N)       , PUBLIC ::   sn_icwmsk
40
41   REAL(wp) :: rn_rztop, rn_rzsta, rn_rzovt
42
43   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_qmelt
44   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_qsubg
45   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_ovt
46   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_ztop
47   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_zovt
48   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf_icwmsk
49
50   INTEGER , PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: nk_top, nk_ovt, nk_bot, nk_sta
51   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: h_neg, h_pos, h_sta
52   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rztop, rzovt, rzsta
53   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rovt 
54   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: qtot, qtot_b 
55   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: ricwmsk
56   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: hdivicw_b, hdivicw 
57   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: ricw_tsc_b, ricw_tsc
58
59   REAL(wp), PUBLIC, SAVE ::   lfusicw= 0.334e6_wp    ! phycst ?
60   REAL(wp), PUBLIC, SAVE ::   E0, K0, r1_E0K0 ! entrainment
61   REAL(wp), PUBLIC, SAVE ::   FPa, FPb, FPc   ! Freezing point
62   REAL(wp), PUBLIC, SAVE ::   GamT, GamTS0, zeropt0
63   REAL(wp), PUBLIC, SAVE ::   ricwpol1, ricwpol2, ricwpol3, ricwpol4 , ricwpol5 , ricwpol6, &
64   &                           ricwpol7, ricwpol8, ricwpol9, ricwpol10, ricwpol11              ! poly coef
65   !! * Substitutions
66#  include "domzgr_substitute.h90"
67#  include "vectopt_loop_substitute.h90"
68   !!----------------------------------------------------------------------
69   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
70   !! $Id: trasbc.F90 4726 2014-07-23 16:27:21Z mathiot $
71   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
72   !!----------------------------------------------------------------------
73CONTAINS
74
75   INTEGER FUNCTION tra_icw_alloc()
76      !!----------------------------------------------------------------------
77      !!                ***  ROUTINE tra_icw_alloc  ***
78      !!----------------------------------------------------------------------
79      ALLOCATE( h_neg  (jpi,jpj), h_pos  (jpi,jpj), h_sta (jpi,jpj), &
80         &      nk_top (jpi,jpj), nk_ovt (jpi,jpj), nk_bot(jpi,jpj), nk_sta(jpi,jpj),                       &
81         &      rovt   (jpi,jpj), qtot   (jpi,jpj), qtot_b(jpi,jpj),                                        & 
82         &      rztop  (jpi,jpj), rzovt  (jpi,jpj), rzsta (jpi,jpj),                                        &
83         &      ricwmsk(jpi,jpj), hdivicw(jpi,jpj,jpk), hdivicw_b(jpi,jpj,jpk),                             &
84         &      ricw_tsc_b(jpi,jpj,jpk,2), ricw_tsc (jpi,jpj,jpk,2), STAT=tra_icw_alloc )
85         !
86      IF( lk_mpp            )   CALL mpp_sum ( tra_icw_alloc )
87      IF( tra_icw_alloc > 0 )   CALL ctl_warn('tra_icw_alloc: allocation of arrays failed')
88   END FUNCTION tra_icw_alloc
89
90
91   SUBROUTINE tra_icw_init
92      !!----------------------------------------------------------------------
93      !!                  ***  ROUTINE tra_icw_init  ***
94      !!
95      !! ** Purpose :   Initialisation of parametrisation data
96      !!
97      !! ** Method  : - read the runoff namtra_icw namelist
98      !!
99      !! ** Action  : - read parameters
100      !!----------------------------------------------------------------------------------
101      CHARACTER(len=32) ::   rn_icw_file   ! runoff file name
102      INTEGER           ::   ios           ! Local integer output status for namelist read
103      INTEGER           ::   ierror,inum   ! temporary integer
104      !
105      NAMELIST/namtra_icw/ cn_dir  , ln_traicw, ln_plumes, ln_ovt ,                        &
106        &                  sn_qmelt, sn_qsubg , sn_ovt   , sn_ztop, sn_zovt, sn_icwmsk, rn_rztop, rn_rzovt, rn_rzsta
107      !!----------------------------------------------------------------------------------
108      !
109      !                                   ! ============
110      !                                   !   Namelist
111      !                                   ! ============
112      !
113      REWIND( numnam_ref )              ! Namelist namtra_icw in reference namelist : Runoffs
114      READ  ( numnam_ref, namtra_icw, IOSTAT = ios, ERR = 901)
115901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_icw in reference namelist', lwp )
116
117      REWIND( numnam_cfg )              ! Namelist namtra_icw in configuration namelist : Runoffs
118      READ  ( numnam_cfg, namtra_icw, IOSTAT = ios, ERR = 902 )
119902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_icw in configuration namelist', lwp )
120      IF(lwm) WRITE ( numond, namtra_icw )
121      !
122      !                                         ! Control print
123      IF(lwp) THEN
124         WRITE(numout,*)
125         WRITE(numout,*) 'tra_icw : parametrisation of melting along calving face '
126         WRITE(numout,*) '~~~~~~~ '
127         WRITE(numout,*) ' use plumes model poly   = ', ln_plumes
128         WRITE(numout,*) ' use prescription of ovt = ', ln_ovt
129      ENDIF
130      !
131      !                                         !==  allocate runoff arrays
132      IF( tra_icw_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_icw_alloc : unable to allocate arrays' )
133      IF (ln_traicw) THEN
134
135        IF( .NOT. ln_plumes ) THEN           !* set up file structure
136           IF(lwp) WRITE(numout,*)
137           IF(lwp) WRITE(numout,*) '        icw parametrisation variables read in a file'
138           ALLOCATE( sf_qmelt(1), STAT=ierror )
139           IF( ierror > 0 ) THEN
140              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_qmelt structure' )   ;   RETURN
141           ENDIF
142           ALLOCATE( sf_qmelt(1)%fnow(jpi,jpj,1)   )
143           IF( sn_qmelt%ln_tint )ALLOCATE( sf_qmelt(1)%fdta(jpi,jpj,1,2) )
144           !                                        ! fill sf_chl with sn_chl and control print
145           CALL fld_fill( sf_qmelt, (/ sn_qmelt /), cn_dir, 'tra_icw_init',   &
146              &                                         'ice wall parametrisation function of read qmelt', 'namtra_icw' )
147               !
148           ALLOCATE( sf_ovt(1), STAT=ierror )
149           IF( ierror > 0 ) THEN
150              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_ovt structure' )   ;   RETURN
151           ENDIF
152           ALLOCATE( sf_ovt(1)%fnow(jpi,jpj,1)   )
153           IF( sn_ovt%ln_tint )ALLOCATE( sf_ovt(1)%fdta(jpi,jpj,1,2) )
154           !                                        ! fill sf_chl with sn_chl and control print
155           CALL fld_fill( sf_ovt, (/ sn_ovt /), cn_dir, 'tra_icw_init',   &
156              &                                         'ice wall parametrisation function of read ovt', 'namtra_icw' )
157
158           ALLOCATE( sf_ztop(1), STAT=ierror )
159           IF( ierror > 0 ) THEN
160              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_ztop structure' )   ;   RETURN
161           ENDIF
162           ALLOCATE( sf_ztop(1)%fnow(jpi,jpj,1)   )
163           IF( sn_ztop%ln_tint )ALLOCATE( sf_ztop(1)%fdta(jpi,jpj,1,2) )
164           !                                        ! fill sf_chl with sn_chl and control print
165           CALL fld_fill( sf_ztop, (/ sn_ztop /), cn_dir, 'tra_icw_init',   &
166              &                                         'ice wall parametrisation function of read ztop', 'namtra_icw' )
167       
168           ALLOCATE( sf_zovt(1), STAT=ierror )
169           IF( ierror > 0 ) THEN
170              CALL ctl_stop( 'tra_icw_init: unable to allocate sf_zovt structure' )   ;   RETURN
171           ENDIF
172           ALLOCATE( sf_zovt(1)%fnow(jpi,jpj,1)   )
173           IF( sn_zovt%ln_tint )ALLOCATE( sf_zovt(1)%fdta(jpi,jpj,1,2) )
174           !                                        ! fill sf_chl with sn_chl and control print
175           CALL fld_fill( sf_zovt, (/ sn_zovt /), cn_dir, 'tra_icw_init',   &
176              &                                         'ice wall parametrisation function of read zovt', 'namtra_icw' )
177        END IF
178
179        ! set up mask
180        ALLOCATE( sf_icwmsk(1), STAT=ierror )
181        IF( ierror > 0 ) THEN
182           CALL ctl_stop( 'tra_icw_init: unable to allocate sf_icwmsk structure' )   ;   RETURN
183        ENDIF
184        ALLOCATE( sf_icwmsk(1)%fnow(jpi,jpj,1)   )
185        IF( sn_icwmsk%ln_tint )ALLOCATE( sf_icwmsk(1)%fdta(jpi,jpj,1,2) )
186                                                ! fill sf_chl with sn_chl and control print
187        CALL fld_fill( sf_icwmsk, (/ sn_icwmsk /), cn_dir, 'tra_icw_init',   &
188        &                                         'ice wall parametrisation function of read icwmsk', 'namtra_icw' )
189         
190        IF ( .NOT. ln_plumes ) THEN
191           ! read depth of the top of the plumes
192           rn_icw_file = TRIM( cn_dir )//TRIM( sn_ztop%clname )
193           CALL iom_open ( rn_icw_file, inum )                       
194           CALL iom_get  ( inum, jpdom_data, sn_ztop%clvar, rztop ) 
195           CALL iom_close( inum )                                   
196           ! read depth of the top of the plumes
197           rn_icw_file = TRIM( cn_dir )//TRIM( sn_zovt%clname )
198           CALL iom_open ( rn_icw_file, inum )                       
199           CALL iom_get  ( inum, jpdom_data, sn_zovt%clvar, rzovt ) 
200           CALL iom_close( inum )                                   
201        ELSE
202           !%melt_param_pierre
203!           E0  =  3.60e-2 ! Entrainment coefficient
204!           K0  =  1.00e-3 !
205!           FPa = -5.73e-2 ! Seawater freezing point slope
206!           FPb =  8.32e-2 ! Seawater freezing point offset
207!           FPc =  7.61e-4 ! Depth dependance of freezing point
208!           GamT   = 3.479e-4 !
209!           GamTS0 = 1.897e-4 !
210!           zeropt0 = 0.56    !
211!           ricwpol1 =  2.3109067e5/10._wp ! poly coef 1
212!           ricwpol2 = -6.2399669e5/ 9._wp ! poly coef 2
213!           ricwpol3 =  7.1371620e5/ 8._wp ! poly coef 3
214!           ricwpol4 = -4.5038638e5/ 7._wp ! poly coef 4
215!           ricwpol5 =  1.7122183e5/ 6._wp ! poly coef 5
216!           ricwpol6 = -0.4025004e5/ 5._wp ! poly coef 6
217!           ricwpol7 =  0.0580962e5/ 4._wp ! poly coef 7
218!           ricwpol8 = -0.0050814e5/ 3._wp ! poly coef 8
219!           ricwpol9 =  0.0002680e5/ 2._wp ! poly coef 9
220!           ricwpol10=  0.0000005e5/ 1._wp ! poly coef 10
221!           ricwpol11=  0.0
222           !%melt_param_pierre_plumes_lab
223!           E0  =  7.20e-2 ! Entrainment coefficient
224!           K0  =  2.50e-3 !
225!           FPa = -5.73e-2 ! Seawater freezing point slope
226!           FPb =  8.32e-2 ! Seawater freezing point offset
227!           FPc =  7.61e-4 ! Depth dependance of freezing point
228!           GamT   = 1.100e-3 !
229!           GamTS0 = 6.000e-4 !
230!           zeropt0 = 0.56    !
231!           ricwpol1 =  1.0412536e6/10._wp ! poly coef 1
232!           ricwpol2 = -2.8114553e6/ 9._wp ! poly coef 2
233!           ricwpol3 =  3.2153640e6/ 8._wp ! poly coef 3
234!           ricwpol4 = -2.0287589e6/ 7._wp ! poly coef 4
235!           ricwpol5 =  0.7711177e6/ 6._wp ! poly coef 5
236!           ricwpol6 = -0.1812180e6/ 5._wp ! poly coef 6
237!           ricwpol7 =  0.0261438e6/ 4._wp ! poly coef 7
238!           ricwpol8 = -0.0022843e6/ 3._wp ! poly coef 8
239!           ricwpol9 =  0.0001203e6/ 2._wp ! poly coef 9
240!           ricwpol10=  0.0000002e6/ 1._wp ! poly coef 10
241!           ricwpol11=  0.0
242           !%melt_param_pierre_plumes_lab
243           E0  =  3.60e-2 ! Entrainment coefficient
244           K0  =  2.50e-3 !
245           FPa = -5.73e-2 ! Seawater freezing point slope
246           FPb =  8.32e-2 ! Seawater freezing point offset
247           FPc =  7.61e-4 ! Depth dependance of freezing point
248           GamT   = 1.100e-3 !
249           GamTS0 = 6.000e-4 !
250           zeropt0 = 0.56    !
251           ricwpol1 =  0.7373182e6/10._wp ! poly coef 1
252           ricwpol2 = -1.9907329e6/ 9._wp ! poly coef 2
253           ricwpol3 =  2.2766835e6/ 8._wp ! poly coef 3
254           ricwpol4 = -1.4364516e6/ 7._wp ! poly coef 4
255           ricwpol5 =  0.5459653e6/ 6._wp ! poly coef 5
256           ricwpol6 = -0.1282985e6/ 5._wp ! poly coef 6
257           ricwpol7 =  0.0185074e6/ 4._wp ! poly coef 7
258           ricwpol8 = -0.0016168e6/ 3._wp ! poly coef 8
259           ricwpol9 =  0.0000851e6/ 2._wp ! poly coef 9
260           ricwpol10=  0.0000001e6/ 1._wp ! poly coef 10
261           ricwpol11=  0.0
262           r1_E0K0 =  1.0/(E0+K0)**(1./3.)
263        END IF
264       
265        ! define subglacial file variable
266        ALLOCATE( sf_qsubg(1), STAT=ierror )
267        IF( ierror > 0 ) THEN
268           CALL ctl_stop( 'tra_icw_init: unable to allocate sf_qsubg structure' )   ;   RETURN
269        ENDIF
270        ALLOCATE( sf_qsubg(1)%fnow(jpi,jpj,1)   )
271        IF( sn_qsubg%ln_tint )ALLOCATE( sf_qsubg(1)%fdta(jpi,jpj,1,2) )
272        !                                        ! fill sf_chl with sn_chl and control print
273        CALL fld_fill( sf_qsubg, (/ sn_qsubg /), cn_dir, 'tra_icw_init',   &
274           &                                         'ice wall parametrisation function of read qmelt', 'namtra_icw' )
275
276        ! read mask
277        rn_icw_file = TRIM( cn_dir )//TRIM( sn_icwmsk%clname )
278        CALL iom_open ( rn_icw_file, inum )                       
279        CALL iom_get  ( inum, jpdom_data, sn_icwmsk%clvar, ricwmsk ) 
280        CALL iom_close( inum )                                   
281
282      END IF
283
284      ! initialisation
285      nk_top(:,:)=jpk     ; nk_ovt(:,:)=jpk     ; nk_bot(:,:)=jpk ; nk_sta(:,:)=jpk
286!      rztop(:,:)=0.0      ; rzovt(:,:)=0.0     ; rzsta(:,:)=0.0
287      qtot (:,:) = 0.0_wp ; rovt (:,:) = 0.0_wp ;  qtot_b(:,:)=0.0
288      ricw_tsc(:,:,:,:) = 0.0_wp ; ricw_tsc_b(:,:,:,:) = 0.0_wp
289      h_pos(:,:) = 1.0_wp ; h_neg(:,:) = 1.0_wp ; h_sta(:,:) = 1.0_wp
290      hdivicw(:,:,:) = 0.0_wp
291
292   END SUBROUTINE tra_icw_init
293
294   SUBROUTINE tra_icw ( kt )
295      !!----------------------------------------------------------------------
296      !!                  ***  ROUTINE tra_icw  ***
297      !!
298      !! ** Purpose :   correct T/S
299      !!
300      !! ** Method  :   apply "negative runoff" at the bottom and positive runoff at the top
301      !!
302      !! ** Action  :   apply the heat and salt flux when we impose the overturning
303      !!                along glacier termini face
304      !!----------------------------------------------------------------------
305      INTEGER, INTENT(in) ::   kt          ! ocean time step
306      !
307      INTEGER  ::   ji, jj, jk             ! dummy loop indices
308      !!----------------------------------------------------------------------
309      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztfrz_melt, ztfrz_sub   ! freezing point used for temperature correction
310      REAL(wp), DIMENSION(:,:)  , POINTER       ::   qmelt , qsubg
311      REAL(wp), DIMENSION(:,:)  , POINTER       ::   zglDe, zglDUf, zTw, zSw, zRw
312      REAL(wp), DIMENSION(:,:)  , POINTER       ::   zdeldfp, zGamTS, zxscale, zmscale, zmfscale
313      REAL(wp), DIMENSION(:,:)  , POINTER       ::   zmeltgl, zfluxgl, zxscalegl
314      REAL(wp), DIMENSION(:,:,:), POINTER       ::   zmeltfluxres, zmf, zPf, zMelt
315      REAL(wp), DIMENSION(:,:,:), POINTER       ::   zricw_tsc_sum
316      REAL(wp), DIMENSION(2)                    ::   zTSw
317      REAL(wp) :: zpress, zmeltfres, zrespts
318
319      !
320      CALL wrk_alloc( jpi,jpj,   ztfrz_melt, ztfrz_sub)
321      CALL wrk_alloc( jpi,jpj,   qmelt , qsubg )
322      CALL wrk_alloc( jpi,jpj,2, zricw_tsc_sum ) 
323      CALL wrk_alloc( jpi,jpj, zglDe, zglDUf, zTw, zSw, zRw)
324      CALL wrk_alloc( jpi,jpj, zdeldfp, zGamTS, zxscale, zmscale, zmfscale )
325      CALL wrk_alloc( jpi,jpj, zmeltgl, zfluxgl, zxscalegl)
326      CALL wrk_alloc( jpi,jpj,jpk, zmeltfluxres, zmf, zPf, zMelt) 
327     
328      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
329         !                                         ! ---------------------------------------- !
330         ricw_tsc_b(:,:,:,:) = ricw_tsc(:,:,:,:)               ! where before fields are set at the end of the routine
331         qtot_b(:,:) = qtot(:,:)               ! where before fields are set at the end of the routine
332         ricw_tsc(:,:,:,:) = 0.0_wp
333      ENDIF
334 
335
336      ! read subglacial
337      CALL fld_read ( kt, 1, sf_qsubg  )   
338      qsubg(:,:) = sf_qsubg(1)%fnow(:,:,1)
339
340      nk_top(:,:) = 1     ; nk_ovt(:,:) = 1 ; nk_bot(:,:) = mbkt(:,:)
341
342      IF ( .NOT. ln_plumes ) THEN
343         ! all is define in a file define in the init process
344         ! read ovt
345         CALL fld_read ( kt, 1, sf_ovt   )   
346         rovt(:,:)  = sf_ovt(1)%fnow(:,:,1)
347         ! read melting along the face
348         CALL fld_read ( kt, 1, sf_qmelt  )   
349         qmelt(:,:) = sf_qmelt(1)%fnow(:,:,1)
350         ! compute nk_pos and nk_ovt to force the overturning
351         WHERE (ricwmsk == 1)
352            rztop= rn_rztop
353            rzsta= rn_rzsta
354            rzovt= rn_rzovt
355         END WHERE
356         DO jj = 1, jpj
357            DO ji = 1, jpi
358               jk = 1
359               DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rztop(ji,jj) ) ;  jk = jk + 1 ;  END DO
360               nk_top(ji,jj) = jk
361            END DO
362         END DO
363         DO jj = 1, jpj
364            DO ji = 1, jpi
365               jk = 1
366               DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rzovt(ji,jj) ) ;  jk = jk + 1 ;  END DO
367               nk_ovt(ji,jj) = MAX(jk,nk_top(ji,jj)+1)   ! at least 1 level beneath the top outflow level (traicw_div issue)
368            END DO
369         END DO
370         DO jj = 1, jpj
371            DO ji = 1, jpi
372               jk = 1
373               DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rzsta(ji,jj) ) ;  jk = jk + 1 ;  END DO
374               nk_sta(ji,jj) = jk
375            END DO
376         END DO
377      ELSE
378      !----------------------------------------------------------------------%
379      ! Define physical constants and parameters (as defined within plume model).
380      !De   (:)   = fsdepw_n(ji,jj,:)                ! depth [m]
381         DO jj=1,jpj
382            DO ji=1,jpi
383               zglDe (ji,jj) = -fsdepw_n(ji,jj,nk_bot(ji,jj)+1)  ! Grounding line [m]
384               zglDUf(ji,jj) = qsubg(ji,jj) / e2t (ji,jj)        ! subglacial runoff [m2/s] e2t should be in the input file
385               zTw   (ji,jj) = tsn(ji,jj,nk_bot(ji,jj),1) 
386               zSw   (ji,jj) = tsn(ji,jj,nk_bot(ji,jj),2) 
387            END DO
388         END DO
389
390         !----------------------------------------------------------------------%
391         ! Define basic melting plume scales.
392         zdeldfp(:,:)  = (zTw(:,:) - FPa*zSw(:,:)-FPb-FPc*zglDe(:,:))/FPc;
393   
394         zGamTS (:,:)  = GamT*(0.545+3.5e-5*zdeldfp(:,:)*E0/(GamTS0+E0))
395         zxscale(:,:)  = zeropt0*(zGamTS(:,:)+E0)/(zeropt0*zGamTS(:,:)+E0)
396         zmscale(:,:)  = 3.16e-7*SQRT(1._wp/(K0+E0))*SQRT(zGamTS(:,:)/(zGamTS(:,:)+E0))    &
397         &              * (E0/(zGamTS(:,:)+E0))*((FPc*zdeldfp(:,:))**2.0)
398         zmfscale(:,:) = FPc*zdeldfp(:,:)/(FPc*zdeldfp(:,:)+84.3)*zGamTS(:,:)/(zGamTS(:,:)+E0)
399   
400         !----------------------------------------------------------------------%
401         ! Define scales for freshwater discharge.
402         zmeltgl  (:,:) = 7.9e-3*r1_E0K0*(zGamTS(:,:)*E0/(zGamTS(:,:)+E0))  &
403         &               * (zglDUf(:,:)**(1./3.))*(FPc*zdeldfp(:,:))/1.8_wp
404         zfluxgl  (:,:) = 0.64*E0*r1_E0K0*(zglDUf(:,:)**(1.0/3.0))*zdeldfp(:,:)/zxscale(:,:)/1.2_wp
405         zxscalegl(:,:) = zglDUf(:,:)/(zmeltgl(:,:)+zmscale(:,:));
406
407         !----------------------------------------------------------------------%
408         DO jj=1,jpj
409            DO ji=1,jpi
410               DO jk=1,nk_bot(ji,jj) !useless, if 3d ovt needed have to uncomment it
411         ! Calculate scaled meltwater flux.
412                  zrespts = (-fsdepw_n(ji,jj,jk)-zglDe(ji,jj))*zxscale(ji,jj)/zdeldfp(ji,jj)
413                  zmeltfres = 1._wp-zrespts
414                  zmeltfluxres(ji,jj,jk) =   ricwpol1 * zrespts ** 10._wp  &
415                  &                        + ricwpol2 * zrespts ** 9._wp   &
416                  &                        + ricwpol3 * zrespts ** 8._wp   &
417                  &                        + ricwpol4 * zrespts ** 7._wp   &
418                  &                        + ricwpol5 * zrespts ** 6._wp   &
419                  &                        + ricwpol6 * zrespts ** 5._wp   &
420                  &                        + ricwpol7 * zrespts ** 4._wp   &
421                  &                        + ricwpol8 * zrespts ** 3._wp   &
422                  &                        + ricwpol9 * zrespts ** 2._wp   &
423                  &                        + ricwpol10* zrespts ** 1._wp   &
424                  &                        + ricwpol11* zrespts ** 0._wp
425   
426         ! Derive physical solution
427                  zmf(ji,jj,jk) = zmeltfluxres(ji,jj,jk)*zmscale(ji,jj)*zdeldfp(ji,jj)/zxscale(ji,jj)
428                  zPf(ji,jj,jk) = zmf(ji,jj,jk)*(1._wp-zmeltfres*zmfscale(ji,jj))/(zmeltfres*zmfscale(ji,jj))                      &
429                  &               + zfluxgl(ji,jj)*zrespts
430                  zMelt(ji,jj,jk) = (  zmeltfluxres(ji,jj,jk)*(zmscale(ji,jj)+zmeltgl(ji,jj))                                      &
431                  &                  + zrespts*zmeltgl(ji,jj)*(zxscalegl(ji,jj)*(zdeldfp(ji,jj)+zglDe(ji,jj))/zdeldfp(ji,jj)**2.)) &
432                  &                 / (1._wp+(zxscalegl(ji,jj)*(zdeldfp(ji,jj)+zglDe(ji,jj))/zdeldfp(ji,jj)**2.0))                 &
433                  &                 * zdeldfp(ji,jj)/zxscale(ji,jj) + zglDUf(ji,jj);
434               END DO
435            END DO
436         END DO
437
438         ! initialise stabilisation level and max ovt level
439         nk_sta(:,:) = nk_bot(:,:)
440         nk_ovt(:,:) = nk_bot(:,:)
441         zRw(:,:) = -99.9_wp
442
443         ! compute nk_ovt
444         DO jj=2,jpj
445            DO ji=2,jpi
446
447               ! Calculate freezing temperature
448               qsubg(ji,jj) = zglDUf(ji,jj) * ricwmsk(ji,jj) * rau0 / e2t(ji,jj)
449               zpress = grav*rau0*fsdept(ji,jj,nk_bot(ji,jj))*1.e-04
450               CALL eos_fzp(0.0_wp, ztfrz_sub(ji,jj), zpress)  ! subglacial runoff is fresh water
451
452   
453               DO WHILE (zRw(ji,jj) .LE. rhd(ji,jj,nk_sta(ji,jj)) .AND. (ricwmsk(ji,jj) == 1) .AND. nk_sta(ji,jj) .GE. 1)
454   
455                  ! update nk_sta
456                  nk_sta(ji,jj) = nk_sta(ji,jj)-1
457
458                  ! compute nk_ovt
459                  jk = nk_sta(ji,jj)
460                  DO WHILE ( jk .LT. mbkt(ji,jj) .AND.                                                                    &
461                  &          fsdept_n(ji,jj,jk) < fsdepw(ji,jj,nk_sta(ji,jj)) + MAX(100._wp,fse3t(ji,jj,nk_sta(ji,jj))))
462                    jk = jk + 1
463                  END DO
464                  nk_ovt(ji,jj) = jk
465                 
466                  ! compute overshoot level
467                  jk = 1
468                  DO WHILE ( jk .LT. mbkt(ji,jj) .AND.                                                                    &
469                  &          fsdept_n(ji,jj,jk) < MAX(1._wp, fsdepw(ji,jj,nk_sta(ji,jj)) - 100._wp))
470                    jk = jk + 1
471                  END DO
472                  nk_top(ji,jj) = jk
473
474                  rovt (ji,jj) = zPf(ji,jj,jk) * rau0 / e2t(ji,jj)
475                  qmelt(ji,jj) = (zMelt(ji,jj,jk) - zglDUf(ji,jj)) * rau0 / e2t(ji,jj)
476                  qtot (ji,jj) =qmelt(ji,jj) + qsubg(ji,jj)
477
478                  ! Calculate freezing temperature
479                  zpress = grav*rau0*fsdept(ji,jj,nk_ovt(ji,jj))*1.e-04
480                  CALL eos_fzp(tsn(ji,jj,jk,jp_sal), ztfrz_melt(ji,jj), zpress)
481
482                  ! compute heat and salt content
483                  zTw(ji,jj) = SUM(tsn(ji,jj,jk:nk_bot(ji,jj),1) * fse3t(ji,jj,jk:nk_bot(ji,jj)))/SUM(fse3t(ji,jj,jk:nk_bot(ji,jj)))
484                  zTw(ji,jj) = ( ( rovt(ji,jj) * zTw(ji,jj) + qmelt(ji,jj) * ztfrz_melt(ji,jj) + qsubg(ji,jj) * ztfrz_sub(ji,jj) ) &
485                     &           - qmelt(ji,jj) * lfusicw * r1_rcp ) / ( rovt(ji,jj) + qtot(ji,jj))
486                  zSw(ji,jj) = SUM(tsn(ji,jj,jk:nk_bot(ji,jj),2) * fse3t(ji,jj,jk:nk_bot(ji,jj)))/SUM(fse3t(ji,jj,jk:nk_bot(ji,jj)))
487                  zSw(ji,jj) = ( rovt(ji,jj) * zSw(ji,jj) ) / ( rovt(ji,jj) + qtot(ji,jj) )
488                  zTSw(1)=zTw(ji,jj)
489                  zTSw(2)=zSw(ji,jj)
490                  CALL eos(zTSw(:), fsdept(ji,jj,nk_sta(ji,jj)), zRw(ji,jj)) ! density at nk_sta to test position
491
492               END DO
493
494               ! compute kn_sta
495               nk_sta(ji,jj) = nk_sta(ji,jj)+1
496
497               ! compute nk_ovt
498               jk = nk_sta(ji,jj)
499               DO WHILE ( jk .LE. mbkt(ji,jj) .AND.                                                                    &
500               &          fsdept_n(ji,jj,jk) < fsdepw(ji,jj,nk_sta(ji,jj)) + MAX(100._wp,fse3t(ji,jj,nk_sta(ji,jj))))
501                 jk = jk + 1
502               END DO
503               nk_ovt(ji,jj) = jk
504
505               ! compute overshoot level
506               jk = 1
507               DO WHILE ( jk .LT. mbkt(ji,jj) .AND.                                                                    &
508               &          fsdept_n(ji,jj,jk) < MAX(1._wp, fsdepw(ji,jj,nk_sta(ji,jj)) - 100._wp))
509                 jk = jk + 1
510               END DO
511               nk_top(ji,jj) = jk
512
513               ! compute melt and ovt
514               qmelt(ji,jj) = (zMelt(ji,jj,nk_top(ji,jj)) - zglDUf(ji,jj)) * e2t(ji,jj)
515               qsubg(ji,jj) =                               zglDUf(ji,jj)  * e2t(ji,jj)
516               rovt (ji,jj) = - zPf(ji,jj,nk_top(ji,jj)) * e2t(ji,jj)
517
518            END DO
519         END DO
520!----------------------------------------------------------------------%
521      END IF
522         
523      ! convert ovt (m3/s) in ovt(kg/m2/s)
524      ! rovt negatif
525      rovt (:,:)=rovt (:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:)
526      qmelt(:,:)=qmelt(:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:)
527      qsubg(:,:)=qsubg(:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:)
528      ! compute total fresh water into the system
529      qtot (:,:) = qmelt(:,:) + qsubg(:,:)
530
531      ! useful to lock the loop in case of no ovt parametrisation
532      IF ( ln_ovt .OR. ln_plumes ) THEN
533      ELSE
534         rovt(:,:)  = 0.0_wp
535         nk_ovt(:,:)=nk_bot(:,:)+1
536      END IF
537
538      ! compute h_pos  and h_neg  to force the overturning
539!      nk_sta = nk_top
540      DO ji=1,jpi
541         DO jj=1,jpj
542            h_pos (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_top(ji,jj):nk_ovt(ji,jj)-1)))
543            h_neg (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_top(ji,jj):nk_bot(ji,jj)  )))
544            h_sta (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_sta(ji,jj):nk_ovt(ji,jj)-1)))
545         END DO
546      END DO
547
548      ! compute heat and salt flux due to the ovt (negative part)
549      ricw_tsc(:,:,:,:) = 0.0_wp
550      DO ji=1,jpi
551         DO jj=1,jpj
552            DO jk=nk_top(ji,jj),nk_bot(ji,jj)
553               ricw_tsc(ji,jj,jk,1) = rovt(ji,jj) * tsn(ji,jj,jk,1) * r1_rau0 / h_neg (ji,jj) * ricwmsk(ji,jj)
554               ricw_tsc(ji,jj,jk,2) = rovt(ji,jj) * tsn(ji,jj,jk,2) * r1_rau0 / h_neg (ji,jj) * ricwmsk(ji,jj)
555            END DO
556         END DO
557      END DO
558
559      ! compute total heat and salt for negative part
560      zricw_tsc_sum(:,:,1) = SUM(ricw_tsc(:,:,:,1) * fse3t_n(:,:,:),3) 
561      zricw_tsc_sum(:,:,2) = SUM(ricw_tsc(:,:,:,2) * fse3t_n(:,:,:),3)
562     
563      ! compute heat and salt flux due to the ovt (positive part)
564      DO jj=1,jpj
565         DO ji=1,jpi
566            ! Calculate freezing temperature
567            zpress = grav*rau0*fsdept(ji,jj,nk_bot(ji,jj))*1.e-04
568            CALL eos_fzp(0.0_wp, ztfrz_sub(ji,jj), zpress)  ! subglacial runoff is fresh water
569
570            DO jk=nk_sta(ji,jj),nk_ovt(ji,jj) - 1
571               ! Calculate freezing temperature
572               zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
573               CALL eos_fzp(tsn(ji,jj,jk,jp_sal), ztfrz_melt(ji,jj), zpress)
574
575               ! compute heat and salt content
576               ricw_tsc(ji,jj,jk,1) = ricw_tsc(ji,jj,jk,1) + ( - zricw_tsc_sum(ji,jj,1)                                    &
577                  &                    + ( qmelt(ji,jj) * ztfrz_melt(ji,jj) + qsubg(ji,jj) * ztfrz_sub(ji,jj) ) * r1_rau0  &
578                  &                    - qmelt(ji,jj) * lfusicw * r1_rau0_rcp                           )                  & 
579                  &                  / h_sta(ji,jj) * ricwmsk(ji,jj)
580               ricw_tsc(ji,jj,jk,2) = ricw_tsc(ji,jj,jk,2) + (- zricw_tsc_sum(ji,jj,2) + qtot(ji,jj) * 0.0_wp * r1_rau0        )  &
581                  &                  / h_sta(ji,jj) * ricwmsk(ji,jj)
582            END DO
583         END DO
584      END DO
585
586      IF(lwp .AND. kt == 1) THEN
587         WRITE(numout,*)
588         WRITE(numout,*) 'tra_icw : parametrisation of melting along calving face '
589         WRITE(numout,*) '~~~~~~~ '
590         WRITE(numout,*) ' rzsta, nk_sta, h_sta        = ', MAXVAL(rzsta*ricwmsk), MAXVAL(nk_sta*ricwmsk), MAXVAL(h_sta*ricwmsk)
591         WRITE(numout,*) ' rzovt, nk_ovt, h_ovt        = ', MAXVAL(rzovt*ricwmsk), MAXVAL(nk_ovt*ricwmsk), MAXVAL(h_neg*ricwmsk)
592      ENDIF
593
594      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
595         !                                             ! ---------------------------------------- !
596     ! need to do something for restart before submission into the trunk
597     !    IF( ln_rstart .AND.    &                               !* Restart: read in restart file
598     !       & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN
599     !       IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields red in the restart file'
600     !       CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b )     ! before runoff
601     !       CALL iom_get( numror, jpdom_autoglo, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) )   ! before heat content of runoff
602     !       CALL iom_get( numror, jpdom_autoglo, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) )   ! before salinity content of runoff
603     !    ELSE                                                   !* no restart: set from nit000 values
604     !       IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000'
605             ricw_tsc_b(:,:,:,:) = ricw_tsc(:,:,:,:)
606             qtot_b(:,:) = qtot(:,:)
607     !    ENDIF
608      ENDIF
609     
610      ! output
611      CALL iom_put('qmelt',  qmelt * r1_rau0 * e1e2t(:,:) )
612      CALL iom_put('qsubg',  qsubg * r1_rau0 * e1e2t(:,:) )
613      CALL iom_put('icwovt', rovt * r1_rau0 * e1e2t(:,:) )
614      CALL iom_put('nktop',  FLOAT(nk_top(:,:)) )
615      CALL iom_put('nksta',  FLOAT(nk_sta(:,:)) )
616      CALL iom_put('nkovt',  FLOAT(nk_ovt(:,:)) )
617
618      ! deallocate array
619      CALL wrk_dealloc( jpi,jpj,   ztfrz_melt, ztfrz_sub )
620      CALL wrk_dealloc( jpi,jpj,   qmelt , qsubg )
621      CALL wrk_dealloc( jpi,jpj,2, zricw_tsc_sum ) 
622      CALL wrk_dealloc( jpi,jpj, zglDe, zglDUf, zTw, zSw, zRw)
623      CALL wrk_dealloc( jpi,jpj, zdeldfp, zGamTS, zxscale, zmscale, zmfscale )
624      CALL wrk_dealloc( jpi,jpj, zmeltgl, zfluxgl, zxscalegl)
625      CALL wrk_dealloc( jpi,jpj,jpk, zmeltfluxres, zmf, zPf, zMelt) 
626   END SUBROUTINE tra_icw
627
628   SUBROUTINE div_icw( phdivn )
629      !!----------------------------------------------------------------------
630      !!                  ***  ROUTINE div_icw  ***
631      !!
632      !! ** Purpose :   update the horizontal divergence with the ovt inflow/outflow
633      !!
634      !! ** Method  :
635      !!                CAUTION : icw contribution is positive (inflow) decreasing the
636      !!                          divergence and expressed in m/s
637      !!
638      !! ** Action  :   phdivn   decreased by the runoff inflow
639      !!----------------------------------------------------------------------
640      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
641      !!
642      INTEGER  ::   ji, jj, jk   ! dummy loop indices
643      !!----------------------------------------------------------------------
644      !
645      hdivicw_b(:,:,:) = hdivicw(:,:,:)
646      hdivicw(:,:,:) = 0.0_wp
647      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
648         DO ji = 1, jpi
649            DO jk = nk_top(ji,jj), nk_bot(ji,jj)                 
650               hdivicw(ji,jj,jk) = hdivicw(ji,jj,jk) -     rovt(ji,jj)                 * r1_rau0 / h_neg(ji,jj)
651            END DO
652            DO jk = nk_sta(ji,jj), nk_ovt(ji,jj)-1                   
653               hdivicw(ji,jj,jk) = hdivicw(ji,jj,jk) - ( - rovt(ji,jj) + qtot(ji,jj) ) * r1_rau0 / h_sta(ji,jj)
654            END DO
655         END DO
656      END DO
657      !
658      phdivn(:,:,:) = phdivn(:,:,:) + 0.5 * ( hdivicw(:,:,:) + hdivicw_b(:,:,:) )
659      !
660   END SUBROUTINE div_icw
661     
662END MODULE traicw
663
Note: See TracBrowser for help on using the repository browser.