MODULE traicw !!============================================================================== !! *** MODULE trasbc *** !! Ocean active tracers: surface boundary condition !!============================================================================== !! History : 02/2015 P. Mathiot : original code !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_icw : parametrisation of ice wall processes !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space domain variables USE in_out_manager ! I/O manager USE fldread ! read input field at current time step USE iom USE wrk_nemo ! Memory Allocation USE timing ! Timing USE eosbn2 IMPLICIT NONE PRIVATE PUBLIC tra_icw ! routine called by step.F90 PUBLIC tra_icw_alloc ! routine call in sbcmod module PUBLIC tra_icw_init ! (PUBLIC for TAM) PUBLIC div_icw ! routine called in sshwzv module LOGICAL, PUBLIC :: ln_plumes LOGICAL, PUBLIC :: ln_ovt LOGICAL, PUBLIC :: ln_traicw CHARACTER(len=100), PUBLIC :: cn_dir !: Root directory for location of icw files TYPE(FLD_N) , PUBLIC :: sn_qmelt TYPE(FLD_N) , PUBLIC :: sn_qsubg TYPE(FLD_N) , PUBLIC :: sn_ovt TYPE(FLD_N) , PUBLIC :: sn_ztop TYPE(FLD_N) , PUBLIC :: sn_zovt TYPE(FLD_N) , PUBLIC :: sn_icwmsk REAL(wp) :: rn_rztop, rn_rzsta, rn_rzovt TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_qmelt TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_qsubg TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_ovt TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_ztop TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_zovt TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf_icwmsk INTEGER , PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: nk_top, nk_ovt, nk_bot, nk_sta REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: h_neg, h_pos, h_sta REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rztop, rzovt, rzsta REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: rovt REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: qtot, qtot_b REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: ricwmsk REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: hdivicw_b, hdivicw REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) :: ricw_tsc_b, ricw_tsc REAL(wp), PUBLIC, SAVE :: lfusicw= 0.334e6_wp ! phycst ? REAL(wp), PUBLIC, SAVE :: E0, K0, r1_E0K0 ! entrainment REAL(wp), PUBLIC, SAVE :: FPa, FPb, FPc ! Freezing point REAL(wp), PUBLIC, SAVE :: GamT, GamTS0, zeropt0 REAL(wp), PUBLIC, SAVE :: ricwpol1, ricwpol2, ricwpol3, ricwpol4 , ricwpol5 , ricwpol6, & & ricwpol7, ricwpol8, ricwpol9, ricwpol10, ricwpol11 ! poly coef !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id: trasbc.F90 4726 2014-07-23 16:27:21Z mathiot $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION tra_icw_alloc() !!---------------------------------------------------------------------- !! *** ROUTINE tra_icw_alloc *** !!---------------------------------------------------------------------- ALLOCATE( h_neg (jpi,jpj), h_pos (jpi,jpj), h_sta (jpi,jpj), & & nk_top (jpi,jpj), nk_ovt (jpi,jpj), nk_bot(jpi,jpj), nk_sta(jpi,jpj), & & rovt (jpi,jpj), qtot (jpi,jpj), qtot_b(jpi,jpj), & & rztop (jpi,jpj), rzovt (jpi,jpj), rzsta (jpi,jpj), & & ricwmsk(jpi,jpj), hdivicw(jpi,jpj,jpk), hdivicw_b(jpi,jpj,jpk), & & ricw_tsc_b(jpi,jpj,jpk,2), ricw_tsc (jpi,jpj,jpk,2), STAT=tra_icw_alloc ) ! IF( lk_mpp ) CALL mpp_sum ( tra_icw_alloc ) IF( tra_icw_alloc > 0 ) CALL ctl_warn('tra_icw_alloc: allocation of arrays failed') END FUNCTION tra_icw_alloc SUBROUTINE tra_icw_init !!---------------------------------------------------------------------- !! *** ROUTINE tra_icw_init *** !! !! ** Purpose : Initialisation of parametrisation data !! !! ** Method : - read the runoff namtra_icw namelist !! !! ** Action : - read parameters !!---------------------------------------------------------------------------------- CHARACTER(len=32) :: rn_icw_file ! runoff file name INTEGER :: ios ! Local integer output status for namelist read INTEGER :: ierror,inum ! temporary integer ! NAMELIST/namtra_icw/ cn_dir , ln_traicw, ln_plumes, ln_ovt , & & sn_qmelt, sn_qsubg , sn_ovt , sn_ztop, sn_zovt, sn_icwmsk, rn_rztop, rn_rzovt, rn_rzsta !!---------------------------------------------------------------------------------- ! ! ! ============ ! ! Namelist ! ! ============ ! REWIND( numnam_ref ) ! Namelist namtra_icw in reference namelist : Runoffs READ ( numnam_ref, namtra_icw, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_icw in reference namelist', lwp ) REWIND( numnam_cfg ) ! Namelist namtra_icw in configuration namelist : Runoffs READ ( numnam_cfg, namtra_icw, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_icw in configuration namelist', lwp ) IF(lwm) WRITE ( numond, namtra_icw ) ! ! ! Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'tra_icw : parametrisation of melting along calving face ' WRITE(numout,*) '~~~~~~~ ' WRITE(numout,*) ' use plumes model poly = ', ln_plumes WRITE(numout,*) ' use prescription of ovt = ', ln_ovt ENDIF ! ! !== allocate runoff arrays IF( tra_icw_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_icw_alloc : unable to allocate arrays' ) IF (ln_traicw) THEN IF( .NOT. ln_plumes ) THEN !* set up file structure IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' icw parametrisation variables read in a file' ALLOCATE( sf_qmelt(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'tra_icw_init: unable to allocate sf_qmelt structure' ) ; RETURN ENDIF ALLOCATE( sf_qmelt(1)%fnow(jpi,jpj,1) ) IF( sn_qmelt%ln_tint )ALLOCATE( sf_qmelt(1)%fdta(jpi,jpj,1,2) ) ! ! fill sf_chl with sn_chl and control print CALL fld_fill( sf_qmelt, (/ sn_qmelt /), cn_dir, 'tra_icw_init', & & 'ice wall parametrisation function of read qmelt', 'namtra_icw' ) ! ALLOCATE( sf_ovt(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'tra_icw_init: unable to allocate sf_ovt structure' ) ; RETURN ENDIF ALLOCATE( sf_ovt(1)%fnow(jpi,jpj,1) ) IF( sn_ovt%ln_tint )ALLOCATE( sf_ovt(1)%fdta(jpi,jpj,1,2) ) ! ! fill sf_chl with sn_chl and control print CALL fld_fill( sf_ovt, (/ sn_ovt /), cn_dir, 'tra_icw_init', & & 'ice wall parametrisation function of read ovt', 'namtra_icw' ) ALLOCATE( sf_ztop(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'tra_icw_init: unable to allocate sf_ztop structure' ) ; RETURN ENDIF ALLOCATE( sf_ztop(1)%fnow(jpi,jpj,1) ) IF( sn_ztop%ln_tint )ALLOCATE( sf_ztop(1)%fdta(jpi,jpj,1,2) ) ! ! fill sf_chl with sn_chl and control print CALL fld_fill( sf_ztop, (/ sn_ztop /), cn_dir, 'tra_icw_init', & & 'ice wall parametrisation function of read ztop', 'namtra_icw' ) ALLOCATE( sf_zovt(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'tra_icw_init: unable to allocate sf_zovt structure' ) ; RETURN ENDIF ALLOCATE( sf_zovt(1)%fnow(jpi,jpj,1) ) IF( sn_zovt%ln_tint )ALLOCATE( sf_zovt(1)%fdta(jpi,jpj,1,2) ) ! ! fill sf_chl with sn_chl and control print CALL fld_fill( sf_zovt, (/ sn_zovt /), cn_dir, 'tra_icw_init', & & 'ice wall parametrisation function of read zovt', 'namtra_icw' ) END IF ! set up mask ALLOCATE( sf_icwmsk(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'tra_icw_init: unable to allocate sf_icwmsk structure' ) ; RETURN ENDIF ALLOCATE( sf_icwmsk(1)%fnow(jpi,jpj,1) ) IF( sn_icwmsk%ln_tint )ALLOCATE( sf_icwmsk(1)%fdta(jpi,jpj,1,2) ) ! fill sf_chl with sn_chl and control print CALL fld_fill( sf_icwmsk, (/ sn_icwmsk /), cn_dir, 'tra_icw_init', & & 'ice wall parametrisation function of read icwmsk', 'namtra_icw' ) IF ( .NOT. ln_plumes ) THEN ! read depth of the top of the plumes rn_icw_file = TRIM( cn_dir )//TRIM( sn_ztop%clname ) CALL iom_open ( rn_icw_file, inum ) CALL iom_get ( inum, jpdom_data, sn_ztop%clvar, rztop ) CALL iom_close( inum ) ! read depth of the top of the plumes rn_icw_file = TRIM( cn_dir )//TRIM( sn_zovt%clname ) CALL iom_open ( rn_icw_file, inum ) CALL iom_get ( inum, jpdom_data, sn_zovt%clvar, rzovt ) CALL iom_close( inum ) ELSE !%melt_param_pierre ! E0 = 3.60e-2 ! Entrainment coefficient ! K0 = 1.00e-3 ! ! FPa = -5.73e-2 ! Seawater freezing point slope ! FPb = 8.32e-2 ! Seawater freezing point offset ! FPc = 7.61e-4 ! Depth dependance of freezing point ! GamT = 3.479e-4 ! ! GamTS0 = 1.897e-4 ! ! zeropt0 = 0.56 ! ! ricwpol1 = 2.3109067e5/10._wp ! poly coef 1 ! ricwpol2 = -6.2399669e5/ 9._wp ! poly coef 2 ! ricwpol3 = 7.1371620e5/ 8._wp ! poly coef 3 ! ricwpol4 = -4.5038638e5/ 7._wp ! poly coef 4 ! ricwpol5 = 1.7122183e5/ 6._wp ! poly coef 5 ! ricwpol6 = -0.4025004e5/ 5._wp ! poly coef 6 ! ricwpol7 = 0.0580962e5/ 4._wp ! poly coef 7 ! ricwpol8 = -0.0050814e5/ 3._wp ! poly coef 8 ! ricwpol9 = 0.0002680e5/ 2._wp ! poly coef 9 ! ricwpol10= 0.0000005e5/ 1._wp ! poly coef 10 ! ricwpol11= 0.0 !%melt_param_pierre_plumes_lab ! E0 = 7.20e-2 ! Entrainment coefficient ! K0 = 2.50e-3 ! ! FPa = -5.73e-2 ! Seawater freezing point slope ! FPb = 8.32e-2 ! Seawater freezing point offset ! FPc = 7.61e-4 ! Depth dependance of freezing point ! GamT = 1.100e-3 ! ! GamTS0 = 6.000e-4 ! ! zeropt0 = 0.56 ! ! ricwpol1 = 1.0412536e6/10._wp ! poly coef 1 ! ricwpol2 = -2.8114553e6/ 9._wp ! poly coef 2 ! ricwpol3 = 3.2153640e6/ 8._wp ! poly coef 3 ! ricwpol4 = -2.0287589e6/ 7._wp ! poly coef 4 ! ricwpol5 = 0.7711177e6/ 6._wp ! poly coef 5 ! ricwpol6 = -0.1812180e6/ 5._wp ! poly coef 6 ! ricwpol7 = 0.0261438e6/ 4._wp ! poly coef 7 ! ricwpol8 = -0.0022843e6/ 3._wp ! poly coef 8 ! ricwpol9 = 0.0001203e6/ 2._wp ! poly coef 9 ! ricwpol10= 0.0000002e6/ 1._wp ! poly coef 10 ! ricwpol11= 0.0 !%melt_param_pierre_plumes_lab E0 = 3.60e-2 ! Entrainment coefficient K0 = 2.50e-3 ! FPa = -5.73e-2 ! Seawater freezing point slope FPb = 8.32e-2 ! Seawater freezing point offset FPc = 7.61e-4 ! Depth dependance of freezing point GamT = 1.100e-3 ! GamTS0 = 6.000e-4 ! zeropt0 = 0.56 ! ricwpol1 = 0.7373182e6/10._wp ! poly coef 1 ricwpol2 = -1.9907329e6/ 9._wp ! poly coef 2 ricwpol3 = 2.2766835e6/ 8._wp ! poly coef 3 ricwpol4 = -1.4364516e6/ 7._wp ! poly coef 4 ricwpol5 = 0.5459653e6/ 6._wp ! poly coef 5 ricwpol6 = -0.1282985e6/ 5._wp ! poly coef 6 ricwpol7 = 0.0185074e6/ 4._wp ! poly coef 7 ricwpol8 = -0.0016168e6/ 3._wp ! poly coef 8 ricwpol9 = 0.0000851e6/ 2._wp ! poly coef 9 ricwpol10= 0.0000001e6/ 1._wp ! poly coef 10 ricwpol11= 0.0 r1_E0K0 = 1.0/(E0+K0)**(1./3.) END IF ! define subglacial file variable ALLOCATE( sf_qsubg(1), STAT=ierror ) IF( ierror > 0 ) THEN CALL ctl_stop( 'tra_icw_init: unable to allocate sf_qsubg structure' ) ; RETURN ENDIF ALLOCATE( sf_qsubg(1)%fnow(jpi,jpj,1) ) IF( sn_qsubg%ln_tint )ALLOCATE( sf_qsubg(1)%fdta(jpi,jpj,1,2) ) ! ! fill sf_chl with sn_chl and control print CALL fld_fill( sf_qsubg, (/ sn_qsubg /), cn_dir, 'tra_icw_init', & & 'ice wall parametrisation function of read qmelt', 'namtra_icw' ) ! read mask rn_icw_file = TRIM( cn_dir )//TRIM( sn_icwmsk%clname ) CALL iom_open ( rn_icw_file, inum ) CALL iom_get ( inum, jpdom_data, sn_icwmsk%clvar, ricwmsk ) CALL iom_close( inum ) END IF ! initialisation nk_top(:,:)=jpk ; nk_ovt(:,:)=jpk ; nk_bot(:,:)=jpk ; nk_sta(:,:)=jpk ! rztop(:,:)=0.0 ; rzovt(:,:)=0.0 ; rzsta(:,:)=0.0 qtot (:,:) = 0.0_wp ; rovt (:,:) = 0.0_wp ; qtot_b(:,:)=0.0 ricw_tsc(:,:,:,:) = 0.0_wp ; ricw_tsc_b(:,:,:,:) = 0.0_wp h_pos(:,:) = 1.0_wp ; h_neg(:,:) = 1.0_wp ; h_sta(:,:) = 1.0_wp hdivicw(:,:,:) = 0.0_wp END SUBROUTINE tra_icw_init SUBROUTINE tra_icw ( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_icw *** !! !! ** Purpose : correct T/S !! !! ** Method : apply "negative runoff" at the bottom and positive runoff at the top !! !! ** Action : apply the heat and salt flux when we impose the overturning !! along glacier termini face !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time step ! INTEGER :: ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:) , POINTER :: ztfrz_melt, ztfrz_sub ! freezing point used for temperature correction REAL(wp), DIMENSION(:,:) , POINTER :: qmelt , qsubg REAL(wp), DIMENSION(:,:) , POINTER :: zglDe, zglDUf, zTw, zSw, zRw REAL(wp), DIMENSION(:,:) , POINTER :: zdeldfp, zGamTS, zxscale, zmscale, zmfscale REAL(wp), DIMENSION(:,:) , POINTER :: zmeltgl, zfluxgl, zxscalegl REAL(wp), DIMENSION(:,:,:), POINTER :: zmeltfluxres, zmf, zPf, zMelt REAL(wp), DIMENSION(:,:,:), POINTER :: zricw_tsc_sum REAL(wp), DIMENSION(2) :: zTSw REAL(wp) :: zpress, zmeltfres, zrespts ! CALL wrk_alloc( jpi,jpj, ztfrz_melt, ztfrz_sub) CALL wrk_alloc( jpi,jpj, qmelt , qsubg ) CALL wrk_alloc( jpi,jpj,2, zricw_tsc_sum ) CALL wrk_alloc( jpi,jpj, zglDe, zglDUf, zTw, zSw, zRw) CALL wrk_alloc( jpi,jpj, zdeldfp, zGamTS, zxscale, zmscale, zmfscale ) CALL wrk_alloc( jpi,jpj, zmeltgl, zfluxgl, zxscalegl) CALL wrk_alloc( jpi,jpj,jpk, zmeltfluxres, zmf, zPf, zMelt) IF( kt /= nit000 ) THEN ! Swap of forcing fields ! ! ! ---------------------------------------- ! ricw_tsc_b(:,:,:,:) = ricw_tsc(:,:,:,:) ! where before fields are set at the end of the routine qtot_b(:,:) = qtot(:,:) ! where before fields are set at the end of the routine ricw_tsc(:,:,:,:) = 0.0_wp ENDIF ! read subglacial CALL fld_read ( kt, 1, sf_qsubg ) qsubg(:,:) = sf_qsubg(1)%fnow(:,:,1) nk_top(:,:) = 1 ; nk_ovt(:,:) = 1 ; nk_bot(:,:) = mbkt(:,:) IF ( .NOT. ln_plumes ) THEN ! all is define in a file define in the init process ! read ovt CALL fld_read ( kt, 1, sf_ovt ) rovt(:,:) = sf_ovt(1)%fnow(:,:,1) ! read melting along the face CALL fld_read ( kt, 1, sf_qmelt ) qmelt(:,:) = sf_qmelt(1)%fnow(:,:,1) ! compute nk_pos and nk_ovt to force the overturning WHERE (ricwmsk == 1) rztop= rn_rztop rzsta= rn_rzsta rzovt= rn_rzovt END WHERE DO jj = 1, jpj DO ji = 1, jpi jk = 1 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rztop(ji,jj) ) ; jk = jk + 1 ; END DO nk_top(ji,jj) = jk END DO END DO DO jj = 1, jpj DO ji = 1, jpi jk = 1 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rzovt(ji,jj) ) ; jk = jk + 1 ; END DO nk_ovt(ji,jj) = MAX(jk,nk_top(ji,jj)+1) ! at least 1 level beneath the top outflow level (traicw_div issue) END DO END DO DO jj = 1, jpj DO ji = 1, jpi jk = 1 DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdept_n(ji,jj,jk+1) < rzsta(ji,jj) ) ; jk = jk + 1 ; END DO nk_sta(ji,jj) = jk END DO END DO ELSE !----------------------------------------------------------------------% ! Define physical constants and parameters (as defined within plume model). !De (:) = fsdepw_n(ji,jj,:) ! depth [m] DO jj=1,jpj DO ji=1,jpi zglDe (ji,jj) = -fsdepw_n(ji,jj,nk_bot(ji,jj)+1) ! Grounding line [m] zglDUf(ji,jj) = qsubg(ji,jj) / e2t (ji,jj) ! subglacial runoff [m2/s] e2t should be in the input file zTw (ji,jj) = tsn(ji,jj,nk_bot(ji,jj),1) zSw (ji,jj) = tsn(ji,jj,nk_bot(ji,jj),2) END DO END DO !----------------------------------------------------------------------% ! Define basic melting plume scales. zdeldfp(:,:) = (zTw(:,:) - FPa*zSw(:,:)-FPb-FPc*zglDe(:,:))/FPc; zGamTS (:,:) = GamT*(0.545+3.5e-5*zdeldfp(:,:)*E0/(GamTS0+E0)) zxscale(:,:) = zeropt0*(zGamTS(:,:)+E0)/(zeropt0*zGamTS(:,:)+E0) zmscale(:,:) = 3.16e-7*SQRT(1._wp/(K0+E0))*SQRT(zGamTS(:,:)/(zGamTS(:,:)+E0)) & & * (E0/(zGamTS(:,:)+E0))*((FPc*zdeldfp(:,:))**2.0) zmfscale(:,:) = FPc*zdeldfp(:,:)/(FPc*zdeldfp(:,:)+84.3)*zGamTS(:,:)/(zGamTS(:,:)+E0) !----------------------------------------------------------------------% ! Define scales for freshwater discharge. zmeltgl (:,:) = 7.9e-3*r1_E0K0*(zGamTS(:,:)*E0/(zGamTS(:,:)+E0)) & & * (zglDUf(:,:)**(1./3.))*(FPc*zdeldfp(:,:))/1.8_wp zfluxgl (:,:) = 0.64*E0*r1_E0K0*(zglDUf(:,:)**(1.0/3.0))*zdeldfp(:,:)/zxscale(:,:)/1.2_wp zxscalegl(:,:) = zglDUf(:,:)/(zmeltgl(:,:)+zmscale(:,:)); !----------------------------------------------------------------------% DO jj=1,jpj DO ji=1,jpi DO jk=1,nk_bot(ji,jj) !useless, if 3d ovt needed have to uncomment it ! Calculate scaled meltwater flux. zrespts = (-fsdepw_n(ji,jj,jk)-zglDe(ji,jj))*zxscale(ji,jj)/zdeldfp(ji,jj) zmeltfres = 1._wp-zrespts zmeltfluxres(ji,jj,jk) = ricwpol1 * zrespts ** 10._wp & & + ricwpol2 * zrespts ** 9._wp & & + ricwpol3 * zrespts ** 8._wp & & + ricwpol4 * zrespts ** 7._wp & & + ricwpol5 * zrespts ** 6._wp & & + ricwpol6 * zrespts ** 5._wp & & + ricwpol7 * zrespts ** 4._wp & & + ricwpol8 * zrespts ** 3._wp & & + ricwpol9 * zrespts ** 2._wp & & + ricwpol10* zrespts ** 1._wp & & + ricwpol11* zrespts ** 0._wp ! Derive physical solution zmf(ji,jj,jk) = zmeltfluxres(ji,jj,jk)*zmscale(ji,jj)*zdeldfp(ji,jj)/zxscale(ji,jj) zPf(ji,jj,jk) = zmf(ji,jj,jk)*(1._wp-zmeltfres*zmfscale(ji,jj))/(zmeltfres*zmfscale(ji,jj)) & & + zfluxgl(ji,jj)*zrespts zMelt(ji,jj,jk) = ( zmeltfluxres(ji,jj,jk)*(zmscale(ji,jj)+zmeltgl(ji,jj)) & & + zrespts*zmeltgl(ji,jj)*(zxscalegl(ji,jj)*(zdeldfp(ji,jj)+zglDe(ji,jj))/zdeldfp(ji,jj)**2.)) & & / (1._wp+(zxscalegl(ji,jj)*(zdeldfp(ji,jj)+zglDe(ji,jj))/zdeldfp(ji,jj)**2.0)) & & * zdeldfp(ji,jj)/zxscale(ji,jj) + zglDUf(ji,jj); END DO END DO END DO ! initialise stabilisation level and max ovt level nk_sta(:,:) = nk_bot(:,:) nk_ovt(:,:) = nk_bot(:,:) zRw(:,:) = -99.9_wp ! compute nk_ovt DO jj=2,jpj DO ji=2,jpi ! Calculate freezing temperature qsubg(ji,jj) = zglDUf(ji,jj) * ricwmsk(ji,jj) * rau0 / e2t(ji,jj) zpress = grav*rau0*fsdept(ji,jj,nk_bot(ji,jj))*1.e-04 CALL eos_fzp(0.0_wp, ztfrz_sub(ji,jj), zpress) ! subglacial runoff is fresh water 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) ! update nk_sta nk_sta(ji,jj) = nk_sta(ji,jj)-1 ! compute nk_ovt jk = nk_sta(ji,jj) DO WHILE ( jk .LT. mbkt(ji,jj) .AND. & & fsdept_n(ji,jj,jk) < fsdepw(ji,jj,nk_sta(ji,jj)) + MAX(100._wp,fse3t(ji,jj,nk_sta(ji,jj)))) jk = jk + 1 END DO nk_ovt(ji,jj) = jk ! compute overshoot level jk = 1 DO WHILE ( jk .LT. mbkt(ji,jj) .AND. & & fsdept_n(ji,jj,jk) < MAX(1._wp, fsdepw(ji,jj,nk_sta(ji,jj)) - 100._wp)) jk = jk + 1 END DO nk_top(ji,jj) = jk rovt (ji,jj) = zPf(ji,jj,jk) * rau0 / e2t(ji,jj) qmelt(ji,jj) = (zMelt(ji,jj,jk) - zglDUf(ji,jj)) * rau0 / e2t(ji,jj) qtot (ji,jj) =qmelt(ji,jj) + qsubg(ji,jj) ! Calculate freezing temperature zpress = grav*rau0*fsdept(ji,jj,nk_ovt(ji,jj))*1.e-04 CALL eos_fzp(tsn(ji,jj,jk,jp_sal), ztfrz_melt(ji,jj), zpress) ! compute heat and salt content 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))) zTw(ji,jj) = ( ( rovt(ji,jj) * zTw(ji,jj) + qmelt(ji,jj) * ztfrz_melt(ji,jj) + qsubg(ji,jj) * ztfrz_sub(ji,jj) ) & & - qmelt(ji,jj) * lfusicw * r1_rcp ) / ( rovt(ji,jj) + qtot(ji,jj)) 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))) zSw(ji,jj) = ( rovt(ji,jj) * zSw(ji,jj) ) / ( rovt(ji,jj) + qtot(ji,jj) ) zTSw(1)=zTw(ji,jj) zTSw(2)=zSw(ji,jj) CALL eos(zTSw(:), fsdept(ji,jj,nk_sta(ji,jj)), zRw(ji,jj)) ! density at nk_sta to test position END DO ! compute kn_sta nk_sta(ji,jj) = nk_sta(ji,jj)+1 ! compute nk_ovt jk = nk_sta(ji,jj) DO WHILE ( jk .LE. mbkt(ji,jj) .AND. & & fsdept_n(ji,jj,jk) < fsdepw(ji,jj,nk_sta(ji,jj)) + MAX(100._wp,fse3t(ji,jj,nk_sta(ji,jj)))) jk = jk + 1 END DO nk_ovt(ji,jj) = jk ! compute overshoot level jk = 1 DO WHILE ( jk .LT. mbkt(ji,jj) .AND. & & fsdept_n(ji,jj,jk) < MAX(1._wp, fsdepw(ji,jj,nk_sta(ji,jj)) - 100._wp)) jk = jk + 1 END DO nk_top(ji,jj) = jk ! compute melt and ovt qmelt(ji,jj) = (zMelt(ji,jj,nk_top(ji,jj)) - zglDUf(ji,jj)) * e2t(ji,jj) qsubg(ji,jj) = zglDUf(ji,jj) * e2t(ji,jj) rovt (ji,jj) = - zPf(ji,jj,nk_top(ji,jj)) * e2t(ji,jj) END DO END DO !----------------------------------------------------------------------% END IF ! convert ovt (m3/s) in ovt(kg/m2/s) ! rovt negatif rovt (:,:)=rovt (:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:) qmelt(:,:)=qmelt(:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:) qsubg(:,:)=qsubg(:,:) * rau0 / e1e2t(:,:) * ricwmsk(:,:) ! compute total fresh water into the system qtot (:,:) = qmelt(:,:) + qsubg(:,:) ! useful to lock the loop in case of no ovt parametrisation IF ( ln_ovt .OR. ln_plumes ) THEN ELSE rovt(:,:) = 0.0_wp nk_ovt(:,:)=nk_bot(:,:)+1 END IF ! compute h_pos and h_neg to force the overturning ! nk_sta = nk_top DO ji=1,jpi DO jj=1,jpj h_pos (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_top(ji,jj):nk_ovt(ji,jj)-1))) h_neg (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_top(ji,jj):nk_bot(ji,jj) ))) h_sta (ji,jj) = MAX(1._wp,SUM(fse3t_n(ji,jj,nk_sta(ji,jj):nk_ovt(ji,jj)-1))) END DO END DO ! compute heat and salt flux due to the ovt (negative part) ricw_tsc(:,:,:,:) = 0.0_wp DO ji=1,jpi DO jj=1,jpj DO jk=nk_top(ji,jj),nk_bot(ji,jj) ricw_tsc(ji,jj,jk,1) = rovt(ji,jj) * tsn(ji,jj,jk,1) * r1_rau0 / h_neg (ji,jj) * ricwmsk(ji,jj) ricw_tsc(ji,jj,jk,2) = rovt(ji,jj) * tsn(ji,jj,jk,2) * r1_rau0 / h_neg (ji,jj) * ricwmsk(ji,jj) END DO END DO END DO ! compute total heat and salt for negative part zricw_tsc_sum(:,:,1) = SUM(ricw_tsc(:,:,:,1) * fse3t_n(:,:,:),3) zricw_tsc_sum(:,:,2) = SUM(ricw_tsc(:,:,:,2) * fse3t_n(:,:,:),3) ! compute heat and salt flux due to the ovt (positive part) DO jj=1,jpj DO ji=1,jpi ! Calculate freezing temperature zpress = grav*rau0*fsdept(ji,jj,nk_bot(ji,jj))*1.e-04 CALL eos_fzp(0.0_wp, ztfrz_sub(ji,jj), zpress) ! subglacial runoff is fresh water DO jk=nk_sta(ji,jj),nk_ovt(ji,jj) - 1 ! Calculate freezing temperature zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 CALL eos_fzp(tsn(ji,jj,jk,jp_sal), ztfrz_melt(ji,jj), zpress) ! compute heat and salt content ricw_tsc(ji,jj,jk,1) = ricw_tsc(ji,jj,jk,1) + ( - zricw_tsc_sum(ji,jj,1) & & + ( qmelt(ji,jj) * ztfrz_melt(ji,jj) + qsubg(ji,jj) * ztfrz_sub(ji,jj) ) * r1_rau0 & & - qmelt(ji,jj) * lfusicw * r1_rau0_rcp ) & & / h_sta(ji,jj) * ricwmsk(ji,jj) 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 ) & & / h_sta(ji,jj) * ricwmsk(ji,jj) END DO END DO END DO IF(lwp .AND. kt == 1) THEN WRITE(numout,*) WRITE(numout,*) 'tra_icw : parametrisation of melting along calving face ' WRITE(numout,*) '~~~~~~~ ' WRITE(numout,*) ' rzsta, nk_sta, h_sta = ', MAXVAL(rzsta*ricwmsk), MAXVAL(nk_sta*ricwmsk), MAXVAL(h_sta*ricwmsk) WRITE(numout,*) ' rzovt, nk_ovt, h_ovt = ', MAXVAL(rzovt*ricwmsk), MAXVAL(nk_ovt*ricwmsk), MAXVAL(h_neg*ricwmsk) ENDIF IF( kt == nit000 ) THEN ! set the forcing field at nit000 - 1 ! ! ! ---------------------------------------- ! ! need to do something for restart before submission into the trunk ! IF( ln_rstart .AND. & !* Restart: read in restart file ! & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN ! IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields red in the restart file' ! CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b ) ! before runoff ! CALL iom_get( numror, jpdom_autoglo, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) ) ! before heat content of runoff ! CALL iom_get( numror, jpdom_autoglo, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) ) ! before salinity content of runoff ! ELSE !* no restart: set from nit000 values ! IF(lwp) WRITE(numout,*) ' nit000-1 runoff forcing fields set to nit000' ricw_tsc_b(:,:,:,:) = ricw_tsc(:,:,:,:) qtot_b(:,:) = qtot(:,:) ! ENDIF ENDIF ! output CALL iom_put('qmelt', qmelt * r1_rau0 * e1e2t(:,:) ) CALL iom_put('qsubg', qsubg * r1_rau0 * e1e2t(:,:) ) CALL iom_put('icwovt', rovt * r1_rau0 * e1e2t(:,:) ) CALL iom_put('nktop', FLOAT(nk_top(:,:)) ) CALL iom_put('nksta', FLOAT(nk_sta(:,:)) ) CALL iom_put('nkovt', FLOAT(nk_ovt(:,:)) ) ! deallocate array CALL wrk_dealloc( jpi,jpj, ztfrz_melt, ztfrz_sub ) CALL wrk_dealloc( jpi,jpj, qmelt , qsubg ) CALL wrk_dealloc( jpi,jpj,2, zricw_tsc_sum ) CALL wrk_dealloc( jpi,jpj, zglDe, zglDUf, zTw, zSw, zRw) CALL wrk_dealloc( jpi,jpj, zdeldfp, zGamTS, zxscale, zmscale, zmfscale ) CALL wrk_dealloc( jpi,jpj, zmeltgl, zfluxgl, zxscalegl) CALL wrk_dealloc( jpi,jpj,jpk, zmeltfluxres, zmf, zPf, zMelt) END SUBROUTINE tra_icw SUBROUTINE div_icw( phdivn ) !!---------------------------------------------------------------------- !! *** ROUTINE div_icw *** !! !! ** Purpose : update the horizontal divergence with the ovt inflow/outflow !! !! ** Method : !! CAUTION : icw contribution is positive (inflow) decreasing the !! divergence and expressed in m/s !! !! ** Action : phdivn decreased by the runoff inflow !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phdivn ! horizontal divergence !! INTEGER :: ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- ! hdivicw_b(:,:,:) = hdivicw(:,:,:) hdivicw(:,:,:) = 0.0_wp DO jj = 1, jpj ! update the depth over which runoffs are distributed DO ji = 1, jpi DO jk = nk_top(ji,jj), nk_bot(ji,jj) hdivicw(ji,jj,jk) = hdivicw(ji,jj,jk) - rovt(ji,jj) * r1_rau0 / h_neg(ji,jj) END DO DO jk = nk_sta(ji,jj), nk_ovt(ji,jj)-1 hdivicw(ji,jj,jk) = hdivicw(ji,jj,jk) - ( - rovt(ji,jj) + qtot(ji,jj) ) * r1_rau0 / h_sta(ji,jj) END DO END DO END DO ! phdivn(:,:,:) = phdivn(:,:,:) + 0.5 * ( hdivicw(:,:,:) + hdivicw_b(:,:,:) ) ! END SUBROUTINE div_icw END MODULE traicw