Changeset 12165 for NEMO/branches/2019/dev_ASINTER-01-05_merged/src
- Timestamp:
- 2019-12-11T09:27:27+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_ASINTER-01-05_merged/src
- Files:
-
- 1 deleted
- 16 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/ice.F90
r11586 r12165 328 328 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sz_i !: ice salinity [PSS] 329 329 330 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip !: melt pond fraction per grid cell area330 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip !: melt pond concentration 331 331 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip !: melt pond volume per grid cell area [m] 332 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_frac !: melt pond volume per ice area333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_ip !: melt pond thickness[m]334 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond fraction332 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_ip_frac !: melt pond fraction (a_ip/a_i) 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: h_ip !: melt pond depth [m] 334 335 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_ip !: total melt pond concentration 336 336 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hm_ip !: mean melt pond depth [m] 337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_ip !: total melt pond volume per unit area[m]337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_ip !: total melt pond volume per gridcell area [m] 338 338 339 339 !!---------------------------------------------------------------------- -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icectl.F90
r11586 r12165 44 44 PUBLIC ice_prt3D 45 45 46 ! thresold values for conservation46 ! thresold rates for conservation 47 47 ! these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk 48 REAL(wp), PARAMETER :: zchk_m = 1.e-5 ! kg/m2/s <=> 1mm of ice per yearspuriously gained/lost49 REAL(wp), PARAMETER :: zchk_s = 1.e-4 ! g/m2/s <=> 1mm of ice per yearspuriously gained/lost (considering s=10g/kg)50 REAL(wp), PARAMETER :: zchk_t = 3. ! W/m2 <=> 1mm of ice per yearspuriously gained/lost (considering Lf=3e5J/kg)48 REAL(wp), PARAMETER :: zchk_m = 2.5e-7 ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost 49 REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg) 50 REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg) 51 51 52 52 !! * Substitutions … … 68 68 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 69 69 !! It prints in ocean.output if there is a violation of conservation at each time-step 70 !! The thresholds (zchk_m, zchk_s, zchk_t) which determine violations are set to 71 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 70 !! The thresholds (zchk_m, zchk_s, zchk_t) determine violations 72 71 !! For salt and heat thresholds, ice is considered to have a salinity of 10 73 72 !! and a heat content of 3e5 J/kg (=latent heat of fusion) … … 133 132 134 133 ! -- advection scheme is conservative? -- ! 135 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 136 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 134 zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather) 135 zetrp = glob_sum( 'icectl', ( diag_trp_ei + diag_trp_es ) * e1e2t ) ! must be close to 0 (only for Prather) 137 136 138 137 ! ice area (+epsi10 to set a threshold > 0 when there is no ice) … … 157 156 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_amax 158 157 ! check if advection scheme is conservative 159 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 160 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 158 ! only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc) 159 ! so the formulation for conservation is different (and not coded) 160 ! it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now 161 !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 162 ! & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 161 163 ENDIF 162 164 ! … … 173 175 !! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true 174 176 !! It prints in ocean.output if there is a violation of conservation at each time-step 175 !! The thresholds (zchk_m, zchk_s, zchk_t) which determine the violation are set to 176 !! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 177 !! The thresholds (zchk_m, zchk_s, zchk_t) determine the violations 177 178 !! For salt and heat thresholds, ice is considered to have a salinity of 10 178 179 !! and a heat content of 3e5 J/kg (=latent heat of fusion) -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_pra.F90
r10425 r12165 16 16 !! adv_pra_rst : read/write Prather field in ice restart file, or initialized to zero 17 17 !!---------------------------------------------------------------------- 18 USE phycst ! physical constant 18 19 USE dom_oce ! ocean domain 19 20 USE ice ! sea-ice variables 20 21 USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call 22 USE icevar ! sea-ice: operations 21 23 ! 22 24 USE in_out_manager ! I/O manager … … 25 27 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 28 USE lbclnk ! lateral boundary conditions (or mpp links) 27 USE prtctl ! Print control28 29 29 30 IMPLICIT NONE … … 36 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness 37 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! lead fraction39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity 40 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxage, syage, sxxage, syyage, sxyage ! ice age 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sxopw, syopw, sxxopw, syyopw, sxyopw ! open water in sea ice42 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content 43 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content … … 81 81 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 82 82 ! 83 INTEGER :: jk, jl, jt ! dummy loop indices 84 INTEGER :: initad ! number of sub-timestep for the advection 85 REAL(wp) :: zcfl , zusnit ! - - 86 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zarea 87 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0opw 88 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0smi, z0oi 89 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ap , z0vp 90 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0es 91 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0ei 83 INTEGER :: ji,jj, jk, jl, jt ! dummy loop indices 84 INTEGER :: icycle ! number of sub-timestep for the advection 85 REAL(wp) :: zdt ! - - 86 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 87 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 88 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 89 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 90 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 91 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp 92 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 93 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei 92 94 !!---------------------------------------------------------------------- 93 95 ! 94 96 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 95 97 ! 96 ALLOCATE( zarea(jpi,jpj) , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) , & 97 & z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) , & 98 & z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) ) 99 ! 100 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 101 zcfl = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 102 zcfl = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 103 CALL mpp_max( 'icedyn_adv_pra', zcfl ) 98 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 99 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 100 ! this should not affect too much the stability 101 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 102 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) 104 103 105 IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp 106 ELSE ; initad = 1 ; zusnit = 1.0_wp 104 ! non-blocking global communication send zcflnow and receive zcflprv 105 CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 ) 106 107 IF( zcflprv(1) > .5 ) THEN ; icycle = 2 108 ELSE ; icycle = 1 107 109 ENDIF 110 zdt = rdt_ice / REAL(icycle) 108 111 109 zarea(:,:) = e1e2t(:,:) 110 !------------------------- 111 ! transported fields 112 !------------------------- 113 z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:) ! Open water area 114 DO jl = 1, jpl 115 z0snw(:,:,jl) = pv_s (:,:, jl) * e1e2t(:,:) ! Snow volume 116 z0ice(:,:,jl) = pv_i (:,:, jl) * e1e2t(:,:) ! Ice volume 117 z0ai (:,:,jl) = pa_i (:,:, jl) * e1e2t(:,:) ! Ice area 118 z0smi(:,:,jl) = psv_i(:,:, jl) * e1e2t(:,:) ! Salt content 119 z0oi (:,:,jl) = poa_i(:,:, jl) * e1e2t(:,:) ! Age content 120 DO jk = 1, nlay_s 121 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 122 END DO 123 DO jk = 1, nlay_i 124 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 125 END DO 126 IF ( ln_pnd_H12 ) THEN 127 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 128 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 112 ! --- transport --- ! 113 zudy(:,:) = pu_ice(:,:) * e2u(:,:) 114 zvdx(:,:) = pv_ice(:,:) * e1v(:,:) 115 116 DO jt = 1, icycle 117 118 ! record at_i before advection (for open water) 119 zati1(:,:) = SUM( pa_i(:,:,:), dim=3 ) 120 121 ! --- transported fields --- ! 122 DO jl = 1, jpl 123 zarea(:,:,jl) = e1e2t(:,:) 124 z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume 125 z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume 126 z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area 127 z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content 128 z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content 129 DO jk = 1, nlay_s 130 z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content 131 END DO 132 DO jk = 1, nlay_i 133 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 134 END DO 135 IF ( ln_pnd_H12 ) THEN 136 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 137 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 138 ENDIF 139 END DO 140 ! 141 ! !--------------------------------------------! 142 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 143 ! !--------------------------------------------! 144 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 145 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 146 CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 147 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 148 CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 149 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 150 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 151 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 152 CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 153 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 154 ! 155 DO jk = 1, nlay_s !--- snow heat content 156 CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 157 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 158 CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 159 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 160 END DO 161 DO jk = 1, nlay_i !--- ice heat content 162 CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 163 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 164 CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 165 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 166 END DO 167 ! 168 IF ( ln_pnd_H12 ) THEN 169 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 170 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 171 CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 172 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 173 ENDIF 174 ! !--------------------------------------------! 175 ELSE !== even ice time step: adv_y then adv_x ==! 176 ! !--------------------------------------------! 177 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume 178 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) 179 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume 180 CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) 181 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity 182 CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) 183 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration 184 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) 185 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age 186 CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) 187 DO jk = 1, nlay_s !--- snow heat content 188 CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 189 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 190 CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), & 191 & sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) ) 192 END DO 193 DO jk = 1, nlay_i !--- ice heat content 194 CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 195 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 196 CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), & 197 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 198 END DO 199 IF ( ln_pnd_H12 ) THEN 200 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 201 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 202 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 203 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 204 ENDIF 205 ! 129 206 ENDIF 207 208 ! --- Recover the properties from their contents --- ! 209 DO jl = 1, jpl 210 pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 211 pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 212 psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 213 poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 214 pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 215 DO jk = 1, nlay_s 216 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 217 END DO 218 DO jk = 1, nlay_i 219 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 220 END DO 221 IF ( ln_pnd_H12 ) THEN 222 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 223 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 224 ENDIF 225 END DO 226 ! 227 ! derive open water from ice concentration 228 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 229 DO jj = 2, jpjm1 230 DO ji = fs_2, fs_jpim1 231 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 232 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 233 END DO 234 END DO 235 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. ) 236 ! 237 ! --- Ensure non-negative fields --- ! 238 ! Remove negative values (conservation is ensured) 239 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 240 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 241 ! 242 ! --- Ensure snow load is not too big --- ! 243 CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 244 ! 130 245 END DO 131 132 ! !--------------------------------------------!133 IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!134 ! !--------------------------------------------!135 DO jt = 1, initad136 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area137 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )138 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), &139 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )140 DO jl = 1, jpl141 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---142 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )143 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), &144 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )145 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---146 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )147 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), &148 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )149 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---150 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )151 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), &152 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )153 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---154 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )155 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), &156 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )157 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---158 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )159 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), &160 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )161 DO jk = 1, nlay_s !--- snow heat contents ---162 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), &163 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )164 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), &165 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )166 END DO167 DO jk = 1, nlay_i !--- ice heat contents ---168 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), &169 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )170 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), &171 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )172 END DO173 IF ( ln_pnd_H12 ) THEN174 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction --175 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) )176 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), &177 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) )178 CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume --179 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) )180 CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), &181 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) )182 ENDIF183 END DO184 END DO185 ! !--------------------------------------------!186 ELSE !== even ice time step: adv_y then adv_x ==!187 ! !--------------------------------------------!188 DO jt = 1, initad189 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area190 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )191 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:), &192 & sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )193 DO jl = 1, jpl194 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---195 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )196 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl), &197 & sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )198 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---199 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )200 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl), &201 & sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )202 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---203 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )204 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl), &205 & sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )206 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---207 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )208 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi (:,:,jl), sxage(:,:,jl), &209 & sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )210 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---211 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )212 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai (:,:,jl), sxa (:,:,jl), &213 & sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )214 DO jk = 1, nlay_s !--- snow heat contents ---215 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), &216 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )217 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl), &218 & sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )219 END DO220 DO jk = 1, nlay_i !--- ice heat contents ---221 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), &222 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )223 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl), &224 & sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )225 END DO226 IF ( ln_pnd_H12 ) THEN227 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), & !--- melt pond fraction ---228 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) )229 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap (:,:,jl), sxap (:,:,jl), &230 & sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl) )231 CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), & !--- melt pond volume ---232 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) )233 CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp (:,:,jl), sxvp (:,:,jl), &234 & sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl) )235 ENDIF236 END DO237 END DO238 ENDIF239 240 !-------------------------------------------241 ! Recover the properties from their contents242 !-------------------------------------------243 pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1)244 DO jl = 1, jpl245 pv_i (:,:, jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)246 pv_s (:,:, jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)247 psv_i(:,:, jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)248 poa_i(:,:, jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)249 pa_i (:,:, jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)250 DO jk = 1, nlay_s251 pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)252 END DO253 DO jk = 1, nlay_i254 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)255 END DO256 IF ( ln_pnd_H12 ) THEN257 pa_ip (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)258 pv_ip (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)259 ENDIF260 END DO261 !262 DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei )263 246 ! 264 247 IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file … … 267 250 268 251 269 SUBROUTINE adv_x( pd f, put , pcrh, psm , ps0 , &252 SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , & 270 253 & psx, psxx, psy , psyy, psxy ) 271 254 !!---------------------------------------------------------------------- … … 275 258 !! variable on x axis 276 259 !!---------------------------------------------------------------------- 277 REAL(wp) , INTENT(in ) :: pdf ! reduction factor forthe time step278 REAL(wp) 279 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s]280 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psm ! area281 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: ps0 ! field to be advected282 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments283 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments260 REAL(wp) , INTENT(in ) :: pdt ! the time step 261 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 262 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s] 263 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 264 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 265 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 266 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 284 267 !! 285 INTEGER :: ji, jj 286 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! local scalars268 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 269 REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars 287 270 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 288 271 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 291 274 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 292 275 !----------------------------------------------------------------------- 293 294 ! Limitation of moments. 295 296 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 297 298 DO jj = 1, jpj 299 DO ji = 1, jpi 300 zslpmax = MAX( 0._wp, ps0(ji,jj) ) 301 zs1max = 1.5 * zslpmax 302 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) ) 303 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 304 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) ) ) 305 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 306 307 ps0 (ji,jj) = zslpmax 308 psx (ji,jj) = zs1new * rswitch 309 psxx(ji,jj) = zs2new * rswitch 310 psy (ji,jj) = psy (ji,jj) * rswitch 311 psyy(ji,jj) = psyy(ji,jj) * rswitch 312 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 313 END DO 276 ! 277 jcat = SIZE( ps0 , 3 ) ! size of input arrays 278 ! 279 DO jl = 1, jcat ! loop on categories 280 ! 281 ! Limitation of moments. 282 DO jj = 2, jpjm1 283 DO ji = 1, jpi 284 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 285 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 286 ! 287 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 288 zs1max = 1.5 * zslpmax 289 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 290 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 291 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 292 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 293 294 ps0 (ji,jj,jl) = zslpmax 295 psx (ji,jj,jl) = zs1new * rswitch 296 psxx(ji,jj,jl) = zs2new * rswitch 297 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 298 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 299 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 300 END DO 301 END DO 302 303 ! Calculate fluxes and moments between boxes i<-->i+1 304 DO jj = 2, jpjm1 ! Flux from i to i+1 WHEN u GT 0 305 DO ji = 1, jpi 306 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 307 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 308 zalfq = zalf * zalf 309 zalf1 = 1.0 - zalf 310 zalf1q = zalf1 * zalf1 311 ! 312 zfm (ji,jj) = zalf * psm (ji,jj,jl) 313 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 314 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 315 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 316 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 317 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 318 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 319 320 ! Readjust moments remaining in the box. 321 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 322 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 323 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 324 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 325 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 326 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 327 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 328 END DO 329 END DO 330 331 DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0. 332 DO ji = 1, fs_jpim1 333 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 334 zalg (ji,jj) = zalf 335 zalfq = zalf * zalf 336 zalf1 = 1.0 - zalf 337 zalg1 (ji,jj) = zalf1 338 zalf1q = zalf1 * zalf1 339 zalg1q(ji,jj) = zalf1q 340 ! 341 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 342 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 343 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 344 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 345 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 346 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 347 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 348 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 349 END DO 350 END DO 351 352 DO jj = 2, jpjm1 ! Readjust moments remaining in the box. 353 DO ji = fs_2, fs_jpim1 354 zbt = zbet(ji-1,jj) 355 zbt1 = 1.0 - zbet(ji-1,jj) 356 ! 357 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 358 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 359 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 360 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 361 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 362 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 363 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 364 END DO 365 END DO 366 367 ! Put the temporary moments into appropriate neighboring boxes. 368 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0. 369 DO ji = fs_2, fs_jpim1 370 zbt = zbet(ji-1,jj) 371 zbt1 = 1.0 - zbet(ji-1,jj) 372 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 373 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 374 zalf1 = 1.0 - zalf 375 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 376 ! 377 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 378 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 379 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 380 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 381 & + zbt1 * psxx(ji,jj,jl) 382 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 383 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 384 & + zbt1 * psxy(ji,jj,jl) 385 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 386 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 387 END DO 388 END DO 389 390 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0. 391 DO ji = fs_2, fs_jpim1 392 zbt = zbet(ji,jj) 393 zbt1 = 1.0 - zbet(ji,jj) 394 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 395 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 396 zalf1 = 1.0 - zalf 397 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 398 ! 399 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 400 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 401 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 402 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 403 & + ( zalf1 - zalf ) * ztemp ) ) 404 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 405 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 406 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 407 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 408 END DO 409 END DO 410 314 411 END DO 315 412 316 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)317 psm (:,:) = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )318 319 ! Calculate fluxes and moments between boxes i<-->i+1320 DO jj = 1, jpj ! Flux from i to i+1 WHEN u GT 0321 DO ji = 1, jpi322 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )323 zalf = MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)324 zalfq = zalf * zalf325 zalf1 = 1.0 - zalf326 zalf1q = zalf1 * zalf1327 !328 zfm (ji,jj) = zalf * psm (ji,jj)329 zf0 (ji,jj) = zalf * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) ) )330 zfx (ji,jj) = zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )331 zfxx(ji,jj) = zalf * psxx(ji,jj) * zalfq332 zfy (ji,jj) = zalf * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )333 zfxy(ji,jj) = zalfq * psxy(ji,jj)334 zfyy(ji,jj) = zalf * psyy(ji,jj)335 336 ! Readjust moments remaining in the box.337 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj)338 ps0 (ji,jj) = ps0 (ji,jj) - zf0(ji,jj)339 psx (ji,jj) = zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )340 psxx(ji,jj) = zalf1 * zalf1q * psxx(ji,jj)341 psy (ji,jj) = psy (ji,jj) - zfy(ji,jj)342 psyy(ji,jj) = psyy(ji,jj) - zfyy(ji,jj)343 psxy(ji,jj) = zalf1q * psxy(ji,jj)344 END DO345 END DO346 347 DO jj = 1, jpjm1 ! Flux from i+1 to i when u LT 0.348 DO ji = 1, fs_jpim1349 zalf = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj)350 zalg (ji,jj) = zalf351 zalfq = zalf * zalf352 zalf1 = 1.0 - zalf353 zalg1 (ji,jj) = zalf1354 zalf1q = zalf1 * zalf1355 zalg1q(ji,jj) = zalf1q356 !357 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj)358 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )359 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )360 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj) * zalfq361 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )362 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj)363 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj)364 END DO365 END DO366 367 DO jj = 2, jpjm1 ! Readjust moments remaining in the box.368 DO ji = fs_2, fs_jpim1369 zbt = zbet(ji-1,jj)370 zbt1 = 1.0 - zbet(ji-1,jj)371 !372 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )373 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )374 psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )375 psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)376 psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )377 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )378 psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)379 END DO380 END DO381 382 ! Put the temporary moments into appropriate neighboring boxes.383 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0.384 DO ji = fs_2, fs_jpim1385 zbt = zbet(ji-1,jj)386 zbt1 = 1.0 - zbet(ji-1,jj)387 psm(ji,jj) = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)388 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj)389 zalf1 = 1.0 - zalf390 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)391 !392 ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)393 psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)394 psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj) &395 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &396 & + zbt1 * psxx(ji,jj)397 psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj) &398 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj) ) ) &399 & + zbt1 * psxy(ji,jj)400 psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)401 psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)402 END DO403 END DO404 405 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0.406 DO ji = fs_2, fs_jpim1407 zbt = zbet(ji,jj)408 zbt1 = 1.0 - zbet(ji,jj)409 psm(ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )410 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj)411 zalf1 = 1.0 - zalf412 ztemp = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)413 !414 ps0(ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )415 psx(ji,jj) = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )416 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj) &417 & + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) ) &418 & + ( zalf1 - zalf ) * ztemp ) )419 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) &420 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) ) )421 psy(ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )422 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )423 END DO424 END DO425 426 413 !-- Lateral boundary conditions 427 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T', 1., ps0 , 'T', 1. & 428 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 429 & , psxx, 'T', 1., psyy, 'T', 1. & 430 & , psxy, 'T', 1. ) 431 432 IF(ln_ctl) THEN 433 CALL prt_ctl(tab2d_1=psm , clinfo1=' adv_x: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ') 434 CALL prt_ctl(tab2d_1=psx , clinfo1=' adv_x: psx :', tab2d_2=psxx, clinfo2=' psxx : ') 435 CALL prt_ctl(tab2d_1=psy , clinfo1=' adv_x: psy :', tab2d_2=psyy, clinfo2=' psyy : ') 436 CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :') 437 ENDIF 414 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 415 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 416 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 438 417 ! 439 418 END SUBROUTINE adv_x 440 419 441 420 442 SUBROUTINE adv_y( pd f, pvt , pcrh, psm , ps0 , &421 SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , & 443 422 & psx, psxx, psy , psyy, psxy ) 444 423 !!--------------------------------------------------------------------- … … 448 427 !! variable on y axis 449 428 !!--------------------------------------------------------------------- 450 REAL(wp) , INTENT(in ) :: pdf ! reduction factor for thetime step451 REAL(wp) 452 REAL(wp), DIMENSION( jpi,jpj), INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s]453 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psm ! area454 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: ps0 ! field to be advected455 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psx , psy ! 1st moments456 REAL(wp), DIMENSION( jpi,jpj), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments429 REAL(wp) , INTENT(in ) :: pdt ! time step 430 REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0) 431 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s] 432 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area 433 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected 434 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments 435 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments 457 436 !! 458 INTEGER :: ji, jj 459 REAL(wp) :: zs1max, z rdt, zslpmax, ztemp! temporary scalars437 INTEGER :: ji, jj, jl, jcat ! dummy loop indices 438 REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars 460 439 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 461 440 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - … … 464 443 REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - - 465 444 !--------------------------------------------------------------------- 466 467 ! Limitation of moments. 468 469 zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 470 471 DO jj = 1, jpj 472 DO ji = 1, jpi 473 zslpmax = MAX( 0._wp, ps0(ji,jj) ) 474 zs1max = 1.5 * zslpmax 475 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 476 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 477 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) ) ) 478 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 479 ! 480 ps0 (ji,jj) = zslpmax 481 psx (ji,jj) = psx (ji,jj) * rswitch 482 psxx(ji,jj) = psxx(ji,jj) * rswitch 483 psy (ji,jj) = zs1new * rswitch 484 psyy(ji,jj) = zs2new * rswitch 485 psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch 486 END DO 445 ! 446 jcat = SIZE( ps0 , 3 ) ! size of input arrays 447 ! 448 DO jl = 1, jcat ! loop on categories 449 ! 450 ! Limitation of moments. 451 DO jj = 1, jpj 452 DO ji = fs_2, fs_jpim1 453 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 454 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 455 ! 456 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 457 zs1max = 1.5 * zslpmax 458 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 459 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 460 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 461 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 462 ! 463 ps0 (ji,jj,jl) = zslpmax 464 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 465 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 466 psy (ji,jj,jl) = zs1new * rswitch 467 psyy(ji,jj,jl) = zs2new * rswitch 468 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 469 END DO 470 END DO 471 472 ! Calculate fluxes and moments between boxes j<-->j+1 473 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 474 DO ji = fs_2, fs_jpim1 475 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 476 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 477 zalfq = zalf * zalf 478 zalf1 = 1.0 - zalf 479 zalf1q = zalf1 * zalf1 480 ! 481 zfm (ji,jj) = zalf * psm(ji,jj,jl) 482 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 483 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 484 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 485 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 486 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 487 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 488 ! 489 ! Readjust moments remaining in the box. 490 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 491 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 492 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 493 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 494 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 495 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 496 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 497 END DO 498 END DO 499 ! 500 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 501 DO ji = fs_2, fs_jpim1 502 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 503 zalg (ji,jj) = zalf 504 zalfq = zalf * zalf 505 zalf1 = 1.0 - zalf 506 zalg1 (ji,jj) = zalf1 507 zalf1q = zalf1 * zalf1 508 zalg1q(ji,jj) = zalf1q 509 ! 510 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 511 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 512 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 513 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 514 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 515 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 516 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 517 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 518 END DO 519 END DO 520 521 ! Readjust moments remaining in the box. 522 DO jj = 2, jpjm1 523 DO ji = fs_2, fs_jpim1 524 zbt = zbet(ji,jj-1) 525 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 526 ! 527 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 528 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 529 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 530 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 531 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 532 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 533 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 534 END DO 535 END DO 536 537 ! Put the temporary moments into appropriate neighboring boxes. 538 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 539 DO ji = fs_2, fs_jpim1 540 zbt = zbet(ji,jj-1) 541 zbt1 = 1.0 - zbet(ji,jj-1) 542 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 543 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 544 zalf1 = 1.0 - zalf 545 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 546 ! 547 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 548 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 549 & + zbt1 * psy(ji,jj,jl) 550 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 551 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 552 & + zbt1 * psyy(ji,jj,jl) 553 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 554 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 555 & + zbt1 * psxy(ji,jj,jl) 556 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 557 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 558 END DO 559 END DO 560 561 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 562 DO ji = fs_2, fs_jpim1 563 zbt = zbet(ji,jj) 564 zbt1 = 1.0 - zbet(ji,jj) 565 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 566 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 567 zalf1 = 1.0 - zalf 568 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 569 ! 570 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 571 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 572 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 573 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 574 & + ( zalf1 - zalf ) * ztemp ) ) 575 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 576 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 577 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 578 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 579 END DO 580 END DO 581 487 582 END DO 488 583 489 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 490 psm(:,:) = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 491 492 ! Calculate fluxes and moments between boxes j<-->j+1 493 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 494 DO ji = 1, jpi 495 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 496 zalf = MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 497 zalfq = zalf * zalf 498 zalf1 = 1.0 - zalf 499 zalf1q = zalf1 * zalf1 500 ! 501 zfm (ji,jj) = zalf * psm(ji,jj) 502 zf0 (ji,jj) = zalf * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj) + (zalf1-zalf) * psyy(ji,jj) ) ) 503 zfy (ji,jj) = zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 504 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj) 505 zfx (ji,jj) = zalf * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 506 zfxy(ji,jj) = zalfq * psxy(ji,jj) 507 zfxx(ji,jj) = zalf * psxx(ji,jj) 508 ! 509 ! Readjust moments remaining in the box. 510 psm (ji,jj) = psm (ji,jj) - zfm(ji,jj) 511 ps0 (ji,jj) = ps0 (ji,jj) - zf0(ji,jj) 512 psy (ji,jj) = zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 513 psyy(ji,jj) = zalf1 * zalf1q * psyy(ji,jj) 514 psx (ji,jj) = psx (ji,jj) - zfx(ji,jj) 515 psxx(ji,jj) = psxx(ji,jj) - zfxx(ji,jj) 516 psxy(ji,jj) = zalf1q * psxy(ji,jj) 584 !-- Lateral boundary conditions 585 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. & 586 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 587 & , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. ) 588 ! 589 END SUBROUTINE adv_y 590 591 592 SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 593 !!------------------------------------------------------------------- 594 !! *** ROUTINE Hsnow *** 595 !! 596 !! ** Purpose : 1- Check snow load after advection 597 !! 2- Correct pond concentration to avoid a_ip > a_i 598 !! 599 !! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface 600 !! then put the snow excess in the ocean 601 !! 602 !! ** Notes : This correction is crucial because of the call to routine icecor afterwards 603 !! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 604 !! make the snow very thick (if concentration decreases drastically) 605 !! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 606 !!------------------------------------------------------------------- 607 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 608 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip 609 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 610 ! 611 INTEGER :: ji, jj, jl ! dummy loop indices 612 REAL(wp) :: z1_dt, zvs_excess, zfra 613 !!------------------------------------------------------------------- 614 ! 615 z1_dt = 1._wp / pdt 616 ! 617 ! -- check snow load -- ! 618 DO jl = 1, jpl 619 DO jj = 1, jpj 620 DO ji = 1, jpi 621 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 622 ! 623 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 624 ! 625 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 626 ! put snow excess in the ocean 627 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 628 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 629 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 630 ! correct snow volume and heat content 631 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 632 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 633 ENDIF 634 ! 635 ENDIF 636 END DO 517 637 END DO 518 638 END DO 519 639 ! 520 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 521 DO ji = 1, jpi 522 zalf = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 523 zalg (ji,jj) = zalf 524 zalfq = zalf * zalf 525 zalf1 = 1.0 - zalf 526 zalg1 (ji,jj) = zalf1 527 zalf1q = zalf1 * zalf1 528 zalg1q(ji,jj) = zalf1q 529 ! 530 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1) 531 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 532 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 533 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1) * zalfq 534 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 535 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 536 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1) 537 END DO 538 END DO 539 540 ! Readjust moments remaining in the box. 541 DO jj = 2, jpj 542 DO ji = 1, jpi 543 zbt = zbet(ji,jj-1) 544 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 545 ! 546 psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 547 ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 548 psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 549 psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj) 550 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 551 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 552 psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 553 END DO 554 END DO 555 556 ! Put the temporary moments into appropriate neighboring boxes. 557 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 558 DO ji = 1, jpi 559 zbt = zbet(ji,jj-1) 560 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 561 psm(ji,jj) = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 562 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj) 563 zalf1 = 1.0 - zalf 564 ztemp = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 565 ! 566 ps0(ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj) 567 psy(ji,jj) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) & 568 & + zbt1 * psy(ji,jj) 569 psyy(ji,jj) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj) & 570 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 571 & + zbt1 * psyy(ji,jj) 572 psxy(ji,jj) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj) & 573 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) ) & 574 & + zbt1 * psxy(ji,jj) 575 psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 576 psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 577 END DO 578 END DO 579 580 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 581 DO ji = 1, jpi 582 zbt = zbet(ji,jj) 583 zbt1 = ( 1.0 - zbet(ji,jj) ) 584 psm(ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 585 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj) 586 zalf1 = 1.0 - zalf 587 ztemp = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj) 588 ps0 (ji,jj) = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 589 psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 590 psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj) & 591 & + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) & 592 & + ( zalf1 - zalf ) * ztemp ) ) 593 psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj) & 594 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 595 psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 596 psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 597 END DO 598 END DO 599 600 !-- Lateral boundary conditions 601 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T', 1., ps0 , 'T', 1. & 602 & , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes 603 & , psxx, 'T', 1., psyy, 'T', 1. & 604 & , psxy, 'T', 1. ) 605 606 IF(ln_ctl) THEN 607 CALL prt_ctl(tab2d_1=psm , clinfo1=' adv_y: psm :', tab2d_2=ps0 , clinfo2=' ps0 : ') 608 CALL prt_ctl(tab2d_1=psx , clinfo1=' adv_y: psx :', tab2d_2=psxx, clinfo2=' psxx : ') 609 CALL prt_ctl(tab2d_1=psy , clinfo1=' adv_y: psy :', tab2d_2=psyy, clinfo2=' psyy : ') 610 CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :') 611 ENDIF 612 ! 613 END SUBROUTINE adv_y 640 !-- correct pond concentration to avoid a_ip > a_i -- ! 641 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 642 ! 643 END SUBROUTINE Hsnow 614 644 615 645 … … 624 654 ! 625 655 ! !* allocate prather fields 626 ALLOCATE( sxopw(jpi,jpj) , syopw(jpi,jpj) , sxxopw(jpi,jpj) , syyopw(jpi,jpj) , sxyopw(jpi,jpj) , & 627 & sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 656 ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , & 628 657 & sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , & 629 658 & sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , & … … 652 681 !! *** ROUTINE adv_pra_rst *** 653 682 !! 654 !! ** Purpose : Read or write RHGfile in restart file683 !! ** Purpose : Read or write file in restart file 655 684 !! 656 685 !! ** Method : use of IOM library … … 671 700 ! !==========================! 672 701 ! 673 IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sx opw' , ldstop = .FALSE. ) ! file exist: id1>0702 IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0 674 703 ELSE ; id1 = 0 ! no restart: id1=0 675 704 ENDIF … … 689 718 CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn ) 690 719 CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn ) 691 ! ! lead fraction720 ! ! ice concentration 692 721 CALL iom_get( numrir, jpdom_autoglo, 'sxa' , sxa ) 693 722 CALL iom_get( numrir, jpdom_autoglo, 'sya' , sya ) … … 707 736 CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage ) 708 737 CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage ) 709 ! ! open water in sea ice710 CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw )711 CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw )712 CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )713 CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )714 CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )715 738 ! ! snow layers heat content 716 739 DO jk = 1, nlay_s … … 752 775 sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness 753 776 sxsn = 0._wp ; sysn = 0._wp ; sxxsn = 0._wp ; syysn = 0._wp ; sxysn = 0._wp ! snow thickness 754 sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! lead fraction777 sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration 755 778 sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity 756 779 sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age 757 sxopw = 0._wp ; syopw = 0._wp ; sxxopw = 0._wp ; syyopw = 0._wp ; sxyopw = 0._wp ! open water in sea ice758 780 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 759 781 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content … … 786 808 CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn ) 787 809 CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn ) 788 ! ! lead fraction810 ! ! ice concentration 789 811 CALL iom_rstput( iter, nitrst, numriw, 'sxa' , sxa ) 790 812 CALL iom_rstput( iter, nitrst, numriw, 'sya' , sya ) … … 804 826 CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage ) 805 827 CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage ) 806 ! ! open water in sea ice807 CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw )808 CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw )809 CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )810 CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )811 CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )812 828 ! ! snow layers heat content 813 829 DO jk = 1, nlay_s -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_umx.F90
r10945 r12165 83 83 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content 84 84 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration 85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction85 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond concentration 86 86 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 87 87 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content … … 319 319 ! 320 320 !== Ice age ==! 321 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 322 zamsk = 1._wp 323 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 324 & poa_i, poa_i ) 325 ENDIF 321 zamsk = 1._wp 322 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, & 323 & poa_i, poa_i ) 326 324 ! 327 325 !== melt ponds ==! 328 326 IF ( ln_pnd_H12 ) THEN 329 ! fraction327 ! concentration 330 328 zamsk = 1._wp 331 329 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, & … … 1529 1527 !! 3- check whether snow load deplets the snow-ice interface below sea level$ 1530 1528 !! and reduce it by sending the excess in the ocean 1531 !! 4- correct pond fraction to avoid a_ip > a_i1529 !! 4- correct pond concentration to avoid a_ip > a_i 1532 1530 !! 1533 1531 !! ** input : Max thickness of the surrounding 9-points … … 1599 1597 END DO 1600 1598 END DO 1601 ! !-- correct pond fraction to avoid a_ip > a_i1599 ! !-- correct pond concentration to avoid a_ip > a_i 1602 1600 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1603 1601 ! -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_rdgrft.F90
r11587 r12165 86 86 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 87 !!------------------------------------------------------------------- 88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij),&89 & apartf(jpij,0:jpl) , hrmin(jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl),&90 & hrmax (jpij,jpl), hi_hrdg(jpij,jpl) , araft (jpij,jpl),&88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 89 & apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 90 & hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , & 91 91 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 92 … … 137 137 REAL(wp) :: zfac ! local scalar 138 138 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 139 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s)140 139 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 141 140 ! … … 175 174 176 175 ! just needed here 177 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i )178 176 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 179 177 ! needed here and in the iteration loop 178 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i) 180 179 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 181 180 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) … … 187 186 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 188 187 ! 189 ! divergence given by the advection scheme 190 ! (which may not be equal to divu as computed from the velocity field) 191 IF ( ln_adv_Pra ) THEN 192 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 193 ELSEIF( ln_adv_UMx ) THEN 194 zdivu_adv(ji) = zdivu(ji) 195 ENDIF 196 ! 197 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough 198 ! ! to give asum = 1.0 after ridging 188 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough 189 ! ! to give asum = 1.0 after ridging 199 190 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 200 opning(ji) = closing_net(ji) + zdivu _adv(ji)191 opning(ji) = closing_net(ji) + zdivu(ji) 201 192 END DO 202 193 ! … … 215 206 ato_i_1d (ipti) = ato_i_1d (ji) 216 207 closing_net(ipti) = closing_net(ji) 217 zdivu _adv (ipti) = zdivu_adv(ji)208 zdivu (ipti) = zdivu (ji) 218 209 opning (ipti) = opning (ji) 219 210 ENDIF … … 259 250 ELSE 260 251 iterate_ridging = 1 261 zdivu _adv(ji) = zfac * r1_rdtice262 closing_net(ji) = MAX( 0._wp, -zdivu _adv(ji) )263 opning (ji) = MAX( 0._wp, zdivu _adv(ji) )252 zdivu (ji) = zfac * r1_rdtice 253 closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 254 opning (ji) = MAX( 0._wp, zdivu(ji) ) 264 255 ENDIF 265 256 END DO … … 309 300 310 301 ! ! Ice thickness needed for rafting 311 WHERE( pa_i(1:npti,:) > epsi 20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)302 WHERE( pa_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 312 303 ELSEWHERE ; zhi(1:npti,:) = 0._wp 313 304 END WHERE … … 328 319 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 329 320 ! 330 WHERE( zasum(1:npti) > epsi 20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)321 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 331 322 ELSEWHERE ; z1_asum(1:npti) = 0._wp 332 323 END WHERE … … 454 445 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 455 446 ! NOTE: 0 < aksum <= 1 456 WHERE( zaksum(1:npti) > epsi 20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)447 WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 457 448 ELSEWHERE ; closing_gross(1:npti) = 0._wp 458 449 END WHERE … … 537 528 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 538 529 539 IF( a_i_2d(ji,jl1) > epsi 20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)530 IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 540 531 ELSE ; z1_ai(ji) = 0._wp 541 532 ENDIF … … 595 586 ! virtual salt flux to keep salinity constant 596 587 IF( nn_icesal /= 2 ) THEN 597 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) 588 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i 598 589 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice & ! put back sss_m into the ocean 599 590 & - s_i_1d(ji) * vsw * rhoi * r1_rdtice ! and get s_i from the ocean -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/iceitd.F90
r11586 r12165 211 211 CALL itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in 212 212 & g0 (1:npti,1), g1 (1:npti,1), hL (1:npti,1), hR (1:npti,1) ) ! out 213 213 ! 214 214 ! Area lost due to melting of thin ice 215 215 DO ji = 1, npti … … 218 218 ! 219 219 zdh0 = h_i_1d(ji) - h_ib_1d(ji) 220 IF( zdh0 < 0.0 ) THEN ! remove area from category 1220 IF( zdh0 < 0.0 ) THEN ! remove area from category 1 221 221 zdh0 = MIN( -zdh0, hi_max(1) ) 222 222 !Integrate g(1) from 0 to dh0 to estimate area melted … … 226 226 zx1 = zetamax 227 227 zx2 = 0.5 * zetamax * zetamax 228 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 228 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed 229 229 zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i 230 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 231 ! of thin ice (zdamax > 0) 230 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting of thin ice (zdamax > 0) 232 231 ! Remove area, conserving volume 233 232 h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) … … 349 348 DO ji = 1, npti 350 349 ! 351 IF( paice(ji) > epsi10 .AND. phice(ji) > 0._wp) THEN350 IF( paice(ji) > epsi10 .AND. phice(ji) > epsi10 ) THEN 352 351 ! 353 352 ! Initialize hL and hR -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/CRS/README.rst
r10279 r12165 2 2 On line biogeochemistry coarsening 3 3 ********************************** 4 5 .. todo:: 6 7 4 8 5 9 .. contents:: … … 63 67 ! 1, MAX of KZ 64 68 ! 2, MIN of KZ 65 ! 3, 10^(MEAN(LOG(KZ)) 66 ! 4, MEDIANE of KZ 69 ! 3, 10^(MEAN(LOG(KZ)) 70 ! 4, MEDIANE of KZ 67 71 ln_crs_wn = .false. ! wn coarsened (T) or computed using horizontal divergence ( F ) 68 72 ! ! … … 73 77 the north-fold lateral boundary condition (ORCA025, ORCA12, ORCA36, ...). 74 78 - ``nn_msh_crs = 1`` will activate the generation of the coarsened grid meshmask. 75 - ``nn_crs_kz`` is the operator to coarsen the vertical mixing coefficient. 79 - ``nn_crs_kz`` is the operator to coarsen the vertical mixing coefficient. 76 80 - ``ln_crs_wn`` 77 81 … … 80 84 - when ``key_vvl`` is not activated, 81 85 82 - coarsened vertical velocities are computed using horizontal divergence (``ln_crs_wn = .false.``) 86 - coarsened vertical velocities are computed using horizontal divergence (``ln_crs_wn = .false.``) 83 87 - or coarsened vertical velocities are computed with an average operator (``ln_crs_wn = .true.``) 84 88 - ``ln_crs_top = .true.``: should be activated to run BCG model in coarsened space; … … 97 101 98 102 In the [attachment:iodef.xml iodef.xml] file, a "nemo" context is defined and 99 some variable defined in [attachment:file_def.xml file_def.xml] are writted on the ocean-dynamic grid. 103 some variable defined in [attachment:file_def.xml file_def.xml] are writted on the ocean-dynamic grid. 100 104 To write variables on the coarsened grid, and in particular the passive tracers, 101 105 a "nemo_crs" context should be defined in [attachment:iodef.xml iodef.xml] and … … 111 115 interpolated `on-the-fly <http://forge.ipsl.jussieu.fr/nemo/wiki/Users/SetupNewConfiguration/Weight-creator>`_. 112 116 Example of namelist for PISCES : 113 117 114 118 .. code-block:: fortran 115 119 … … 134 138 rn_trfac(14) = 1.0e-06 ! - - - - 135 139 rn_trfac(23) = 7.6e-06 ! - - - - 136 140 137 141 cn_dir = './' ! root directory for the location of the data files 138 142 -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DYN/dynnxt.F90
r10425 r12165 175 175 IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap 176 176 DO jk = 1, jpkm1 177 ub(:,:,jk) = un(:,:,jk) ! ub <-- un 178 vb(:,:,jk) = vn(:,:,jk) 177 179 un(:,:,jk) = ua(:,:,jk) ! un <-- ua 178 180 vn(:,:,jk) = va(:,:,jk) -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flodom.F90
r11413 r12165 433 433 IF( ABS(dlx) > 1.0_wp ) dlx = 1.0_wp 434 434 ! 435 dld = ATAN( DSQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls435 dld = ATAN(SQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls 436 436 flo_dstnce = dld * 1000._wp 437 437 ! -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flowri.F90
r11413 r12165 221 221 clname=TRIM(clname)//".nc" 222 222 223 CALL fliocrfd( clname , (/ 'ntraj' , 't' /), (/ jpnfl , -1/) , numflo )223 CALL fliocrfd( clname , (/'ntraj' , ' t' /), (/ jpnfl , -1/) , numflo ) 224 224 225 225 CALL fliodefv( numflo, 'traj_lon' , (/1,2/), v_t=flio_r8, long_name="Longitude" , units="degrees_east" ) -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LBC/mppini.F90
r11586 r12165 538 538 9401 FORMAT(' ' ,20(' ',i3,' ') ) 539 539 9402 FORMAT(' ',i3,' * ',20(i3,' x',i3,' * ') ) 540 9404 FORMAT(' * ' ,20(' ',i3,' * ') )540 9404 FORMAT(' * ' ,20(' ' ,i4,' * ') ) 541 541 ENDIF 542 542 -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LDF/ldfdyn.F90
r11348 r12165 315 315 DO jj = 1, jpj ! Set local gridscale values 316 316 DO ji = 1, jpi 317 esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2318 esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2317 esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 318 esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 319 319 END DO 320 320 END DO -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OFF/nemogcm.F90
r11348 r12165 114 114 #else 115 115 CALL dta_dyn ( istp ) ! Interpolation of the dynamical fields 116 #endif 117 CALL trc_stp ( istp ) ! time-stepping 118 #if ! defined key_sed_off 116 119 IF( .NOT.ln_linssh ) CALL dta_dyn_swp( istp ) ! swap of sea surface height and vertical scale factors 117 120 #endif 118 CALL trc_stp ( istp ) ! time-stepping119 121 CALL stp_ctl ( istp, indic ) ! Time loop: control and print 120 122 istp = istp + 1 -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/README.rst
r10549 r12165 3 3 *************** 4 4 5 .. todo:: 6 7 8 5 9 .. contents:: 6 :local: 7 8 TOP (Tracers in the Ocean Paradigm) is the NEMO hardwired interface toward biogeochemical models and 9 provide the physical constraints/boundaries for oceanic tracers. 10 It consists of a modular framework to handle multiple ocean tracers, including also a variety of built-in modules. 10 :local: 11 12 TOP (Tracers in the Ocean Paradigm) is the NEMO hardwired interface toward 13 biogeochemical models and provide the physical constraints/boundaries for oceanic tracers. 14 It consists of a modular framework to handle multiple ocean tracers, 15 including also a variety of built-in modules. 11 16 12 17 This component of the NEMO framework allows one to exploit available modules (see below) and 13 18 further develop a range of applications, spanning from the implementation of a dye passive tracer to 14 19 evaluate dispersion processes (by means of MY_TRC), track water masses age (AGE module), 15 assess the ocean interior penetration of persistent chemical compounds (e.g., gases like CFC or even PCBs), 16 up to the full set of equations involving marine biogeochemical cycles. 20 assess the ocean interior penetration of persistent chemical compounds 21 (e.g., gases like CFC or even PCBs), up to the full set of equations involving 22 marine biogeochemical cycles. 17 23 18 24 Structure 19 25 ========= 20 26 21 TOP interface has the following location in the source code ``./src/MBG/`` and27 TOP interface has the following location in the source code :file:`./src/TOP` and 22 28 the following modules are available: 23 29 24 ``TRP`` 25 Interface to NEMO physical core for computing tracers transport 26 27 ``CFC`` 28 Inert carbon tracers (CFC11,CFC12,SF6) 29 30 ``C14`` 31 Radiocarbon passive tracer 32 33 ``AGE`` 34 Water age tracking 35 36 ``MY_TRC`` 37 Template for creation of new modules and external BGC models coupling 38 39 ``PISCES`` 40 Built in BGC model. 41 See [https://www.geosci-model-dev.net/8/2465/2015/gmd-8-2465-2015-discussion.html Aumont et al. (2015)] for 42 a throughout description. | 43 44 The usage of TOP is activated i) by including in the configuration definition the component ``MBG`` and 45 ii) by adding the macro ``key_top`` in the configuration CPP file 46 (see for more details [http://forge.ipsl.jussieu.fr/nemo/wiki/Users "Learn more about the model"]). 30 :file:`TRP` 31 Interface to NEMO physical core for computing tracers transport 32 33 :file:`CFC` 34 Inert carbon tracers (CFC11,CFC12,SF6) 35 36 :file:`C14` 37 Radiocarbon passive tracer 38 39 :file:`AGE` 40 Water age tracking 41 42 :file:`MY_TRC` 43 Template for creation of new modules and external BGC models coupling 44 45 :file:`PISCES` 46 Built in BGC model. See :cite:`gmd-8-2465-2015` for a throughout description. 47 48 The usage of TOP is activated 49 *i)* by including in the configuration definition the component ``TOP`` and 50 *ii)* by adding the macro ``key_top`` in the configuration CPP file 51 (see for more details :forge:`"Learn more about the model" <wiki/Users>`). 47 52 48 53 As an example, the user can refer to already available configurations in the code, … … 51 56 (see also Section 4) . 52 57 53 Note that, since version 4.0, TOP interface core functionalities are activated by means of logical keys and 58 Note that, since version 4.0, 59 TOP interface core functionalities are activated by means of logical keys and 54 60 all submodules preprocessing macros from previous versions were removed. 55 61 … … 57 63 58 64 ``key_iomput`` 59 65 use XIOS I/O 60 66 61 67 ``key_agrif`` 62 68 enable AGRIF coupling 63 69 64 70 ``key_trdtrc`` & ``key_trdmxl_trc`` 65 71 trend computation for tracers 66 72 67 73 Synthetic Workflow 68 74 ================== 69 75 70 A synthetic description of the TOP interface workflow is given below to summarize the steps involved in 71 the computation of biogeochemical and physical trends and their time integration and outputs, 76 A synthetic description of the TOP interface workflow is given below to 77 summarize the steps involved in the computation of biogeochemical and physical trends and 78 their time integration and outputs, 72 79 by reporting also the principal Fortran subroutine herein involved. 73 80 74 **Model initialization (OPA_SRC/nemogcm.F90)** 75 76 call to trc_init (trcini.F90) 77 78 ↳ call trc_nam (trcnam.F90) to initialize TOP tracers and run setting 79 80 ↳ call trc_ini_sms, to initialize each submodule 81 82 ↳ call trc_ini_trp, to initialize transport for tracers 83 84 ↳ call trc_ice_ini, to initialize tracers in seaice 85 86 ↳ call trc_ini_state, read passive tracers from a restart or input data 87 88 ↳ call trc_sub_ini, setup substepping if {{{nn_dttrc /= 1}}} 89 90 **Time marching procedure (OPA_SRC/stp.F90)** 91 92 call to trc_stp.F90 (trcstp.F90) 93 94 ↳ call trc_sub_stp, averaging physical variables for sub-stepping 95 96 ↳ call trc_wri, call XIOS for output of data 97 98 ↳ call trc_sms, compute BGC trends for each submodule 99 100 ↳ call trc_sms_my_trc, includes also surface and coastal BCs trends 101 102 ↳ call trc_trp (TRP/trctrp.F90), compute physical trends 103 104 ↳ call trc_sbc, get trend due to surface concentration/dilution 105 106 ↳ call trc_adv, compute tracers advection 107 108 ↳ call to trc_ldf, compute tracers lateral diffusion 109 110 ↳ call to trc_zdf, vertical mixing and after tracer fields 111 112 ↳ call to trc_nxt, tracer fields at next time step. Lateral Boundary Conditions are solved in here. 113 114 ↳ call to trc_rad, Correct artificial negative concentrations 115 116 ↳ call trc_rst_wri, output tracers restart files 81 Model initialization (:file:`./src/OCE/nemogcm.F90`) 82 ---------------------------------------------------- 83 84 Call to ``trc_init`` subroutine (:file:`./src/TOP/trcini.F90`) to initialize TOP. 85 86 .. literalinclude:: ../../../src/TOP/trcini.F90 87 :language: fortran 88 :lines: 41-86 89 :emphasize-lines: 21,30-32,38-40 90 :caption: ``trc_init`` subroutine 91 92 Time marching procedure (:file:`./src/OCE/step.F90`) 93 ---------------------------------------------------- 94 95 Call to ``trc_stp`` subroutine (:file:`./src/TOP/trcstp.F90`) to compute/update passive tracers. 96 97 .. literalinclude:: ../../../src/TOP/trcstp.F90 98 :language: fortran 99 :lines: 46-125 100 :emphasize-lines: 42,55-57 101 :caption: ``trc_stp`` subroutine 102 103 BGC trends computation for each submodule (:file:`./src/TOP/trcsms.F90`) 104 ------------------------------------------------------------------------ 105 106 .. literalinclude:: ../../../src/TOP/trcsms.F90 107 :language: fortran 108 :lines: 21 109 :caption: :file:`trcsms` snippet 110 111 Physical trends computation (:file:`./src/TOP/TRP/trctrp.F90`) 112 -------------------------------------------------------------- 113 114 .. literalinclude:: ../../../src/TOP/TRP/trctrp.F90 115 :language: fortran 116 :lines: 46-95 117 :emphasize-lines: 17,21,29,33-35 118 :caption: ``trc_trp`` subroutine 117 119 118 120 Namelists walkthrough 119 121 ===================== 120 122 121 namelist_top 122 ------------ 123 124 Here below are listed the features/options of the TOP interface accessible through the namelist_top_ref and 125 modifiable by means of namelist_top_cfg (as for NEMO physical ones). 126 127 Note that ## is used to refer to a number in an array field. 123 :file:`namelist_top` 124 -------------------- 125 126 Here below are listed the features/options of the TOP interface accessible through 127 the :file:`namelist_top_ref` and modifiable by means of :file:`namelist_top_cfg` 128 (as for NEMO physical ones). 129 130 Note that ``##`` is used to refer to a number in an array field. 128 131 129 132 .. literalinclude:: ../../namelists/namtrc_run 133 :language: fortran 130 134 131 135 .. literalinclude:: ../../namelists/namtrc 136 :language: fortran 132 137 133 138 .. literalinclude:: ../../namelists/namtrc_dta 139 :language: fortran 134 140 135 141 .. literalinclude:: ../../namelists/namtrc_adv 142 :language: fortran 136 143 137 144 .. literalinclude:: ../../namelists/namtrc_ldf 145 :language: fortran 138 146 139 147 .. literalinclude:: ../../namelists/namtrc_rad 148 :language: fortran 140 149 141 150 .. literalinclude:: ../../namelists/namtrc_snk 151 :language: fortran 142 152 143 153 .. literalinclude:: ../../namelists/namtrc_dmp 154 :language: fortran 144 155 145 156 .. literalinclude:: ../../namelists/namtrc_ice 157 :language: fortran 146 158 147 159 .. literalinclude:: ../../namelists/namtrc_trd 160 :language: fortran 148 161 149 162 .. literalinclude:: ../../namelists/namtrc_bc 163 :language: fortran 150 164 151 165 .. literalinclude:: ../../namelists/namtrc_bdy 166 :language: fortran 152 167 153 168 .. literalinclude:: ../../namelists/namage 154 155 Two main types of data structure are used within TOP interface to initialize tracer properties (1) and 169 :language: fortran 170 171 Two main types of data structure are used within TOP interface 172 to initialize tracer properties (1) and 156 173 to provide related initial and boundary conditions (2). 157 174 158 **1. TOP tracers initialization**: sn_tracer (namtrc) 175 1. TOP tracers initialization: ``sn_tracer`` (``&namtrc``) 176 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 159 177 160 178 Beside providing name and metadata for tracers, 161 here are also defined the use of initial ({{{sn_tracer%llinit}}}) and 162 boundary ({{{sn_tracer%llsbc, sn_tracer%llcbc, sn_tracer%llobc}}}) conditions. 163 164 In the following, an example of the full structure definition is given for two idealized tracers both with 165 initial conditions given, while the first has only surface boundary forcing and 179 here are also defined the use of initial (``sn_tracer%llinit``) and 180 boundary (``sn_tracer%llsbc, sn_tracer%llcbc, sn_tracer%llobc``) conditions. 181 182 In the following, an example of the full structure definition is given for 183 two idealized tracers both with initial conditions given, 184 while the first has only surface boundary forcing and 166 185 the second both surface and coastal forcings: 167 186 168 187 .. code-block:: fortran 169 188 170 171 172 189 ! ! name ! title of the field ! units ! initial data ! sbc ! cbc ! obc ! 190 sn_tracer(1) = 'TRC1' , 'Tracer 1 Concentration ', ' - ' , .true. , .true., .false., .true. 191 sn_tracer(2) = 'TRC2 ' , 'Tracer 2 Concentration ', ' - ' , .true. , .true., .true. , .false. 173 192 174 193 As tracers in BGC models are increasingly growing, … … 177 196 .. code-block:: fortran 178 197 179 180 181 182 183 184 185 186 198 ! ! name ! title of the field ! units ! initial data ! 199 sn_tracer(1) = 'TRC1' , 'Tracer 1 Concentration ', ' - ' , .true. 200 sn_tracer(2) = 'TRC2 ' , 'Tracer 2 Concentration ', ' - ' , .true. 201 ! sbc 202 sn_tracer(1)%llsbc = .true. 203 sn_tracer(2)%llsbc = .true. 204 ! cbc 205 sn_tracer(2)%llcbc = .true. 187 206 188 207 The data structure is internally initialized by code with dummy names and 189 all initialization/forcing logical fields set to .false. . 190 191 **2. Structures to read input initial and boundary conditions**: namtrc_dta (sn_trcdta), namtrc_bc (sn_trcsbc/sn_trccbc/sn_trcobc) 208 all initialization/forcing logical fields set to ``.false.`` . 209 210 2. Structures to read input initial and boundary conditions: ``&namtrc_dta`` (``sn_trcdta``), ``&namtrc_bc`` (``sn_trcsbc`` / ``sn_trccbc`` / ``sn_trcobc``) 211 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 192 212 193 213 The overall data structure (Fortran type) is based on the general one defined for NEMO core in the SBC component 194 (see details in User Manual SBC Chapteron Input Data specification).195 196 Input fields are prescribed within namtrc_dta (with sn_trcdtastructure),197 while Boundary Conditions are applied to the model by means of namtrc_bc,198 with dedicated structure fields for surface ( sn_trcsbc), riverine (sn_trccbc), and199 lateral open ( sn_trcobc) boundaries.214 (see details in ``SBC`` Chapter of :doc:`Reference Manual <cite>` on Input Data specification). 215 216 Input fields are prescribed within ``&namtrc_dta`` (with ``sn_trcdta`` structure), 217 while Boundary Conditions are applied to the model by means of ``&namtrc_bc``, 218 with dedicated structure fields for surface (``sn_trcsbc``), riverine (``sn_trccbc``), and 219 lateral open (``sn_trcobc``) boundaries. 200 220 201 221 The following example illustrates the data structure in the case of initial condition for 202 a single tracer contained in the file named tracer_1_data.nc (.nc is implicitly assumed in namelist filename), 203 with a doubled initial value, and located in the usr/work/model/inputdata/ folder: 222 a single tracer contained in the file named :file:`tracer_1_data.nc` 223 (``.nc`` is implicitly assumed in namelist filename), 224 with a doubled initial value, and located in the :file:`usr/work/model/inputdata` folder: 204 225 205 226 .. code-block:: fortran 206 227 207 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 208 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 209 sn_trcdta(1) = 'tracer_1_data' , -12 , 'TRC1' , .false. , .true. , 'yearly' , '' , '' , '' 210 rf_trfac(1) = 2.0 211 cn_dir = “usr/work/model/inputdata/” 212 213 Note that, the Lateral Open Boundaries conditions are applied on the segments defined for the physical core of NEMO 214 (see BDY description in the User Manual). 215 216 namelist_trc 217 ------------ 218 219 Here below the description of namelist_trc_ref used to handle Carbon tracers modules, namely CFC and C14. 220 221 |||| &'''namcfc''' ! CFC || 222 223 |||| &'''namc14_typ''' ! C14 - type of C14 tracer, default values of C14/C and pco2 || 224 225 |||| &'''namc14_sbc''' ! C14 - surface BC || 226 227 |||| &'''namc14_fcg''' ! files & dates || 228 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 229 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 230 sn_trcdta(1) = 'tracer_1_data' , -12 , 'TRC1' , .false. , .true. , 'yearly' , '' , '' , '' 231 rf_trfac(1) = 2.0 232 cn_dir = 'usr/work/model/inputdata/' 233 234 Note that, the Lateral Open Boundaries conditions are applied on 235 the segments defined for the physical core of NEMO 236 (see ``BDY`` description in the :doc:`Reference Manual <cite>`). 237 238 :file:`namelist_trc` 239 -------------------- 240 241 Here below the description of :file:`namelist_trc_ref` used to handle Carbon tracers modules, 242 namely CFC and C14. 243 244 .. literalinclude:: ../../../cfgs/SHARED/namelist_trc_ref 245 :language: fortran 246 :lines: 7,17,26,34 247 :caption: :file:`namelist_trc_ref` snippet 228 248 229 249 ``MY_TRC`` interface for coupling external BGC models 230 250 ===================================================== 231 251 232 The generalized interface is pivoted on MY_TRC module that contains template files to build the coupling between 252 The generalized interface is pivoted on MY_TRC module that contains template files to 253 build the coupling between 233 254 NEMO and any external BGC model. 234 255 235 The call to MY_TRC is activated by setting ``ln_my_trc = .true.`` (in namtrc)256 The call to MY_TRC is activated by setting ``ln_my_trc = .true.`` (in ``&namtrc``) 236 257 237 258 The following 6 fortran files are available in MY_TRC with the specific purposes here described. 238 259 239 ``par_my_trc.F90`` 240 This module allows to define additional arrays and public variables to be used within the MY_TRC interface 241 242 ``trcini_my_trc.F90`` 243 Here are initialized user defined namelists and the call to the external BGC model initialization procedures to 244 populate general tracer array (trn and trb). Here are also likely to be defined suport arrays related to 245 system metrics that could be needed by the BGC model. 246 247 ``trcnam_my_trc.F90`` 248 This routine is called at the beginning of trcini_my_trc and should contain the initialization of 249 additional namelists for the BGC model or user-defined code. 250 251 ``trcsms_my_trc.F90`` 252 The routine performs the call to Boundary Conditions and its main purpose is to 253 contain the Source-Minus-Sinks terms due to the biogeochemical processes of the external model. 254 Be aware that lateral boundary conditions are applied in trcnxt routine. 255 IMPORTANT: the routines to compute the light penetration along the water column and 256 the tracer vertical sinking should be defined/called in here, as generalized modules are still missing in 257 the code. 258 259 ``trcice_my_trc.F90`` 260 Here it is possible to prescribe the tracers concentrations in the seaice that will be used as 261 boundary conditions when ice melting occurs (nn_ice_tr =1 in namtrc_ice). 262 See e.g. the correspondent PISCES subroutine. 263 264 ``trcwri_my_trc.F90`` 265 This routine performs the output of the model tracers (only those defined in namtrc) using IOM module 266 (see Manual Chapter “Output and Diagnostics”). 267 It is possible to place here the output of additional variables produced by the model, 268 if not done elsewhere in the code, using the call to iom_put. 260 :file:`par_my_trc.F90` 261 This module allows to define additional arrays and public variables to 262 be used within the MY_TRC interface 263 264 :file:`trcini_my_trc.F90` 265 Here are initialized user defined namelists and 266 the call to the external BGC model initialization procedures to populate general tracer array 267 (``trn`` and ``trb``). 268 Here are also likely to be defined support arrays related to system metrics that 269 could be needed by the BGC model. 270 271 :file:`trcnam_my_trc.F90` 272 This routine is called at the beginning of ``trcini_my_trc`` and 273 should contain the initialization of additional namelists for the BGC model or user-defined code. 274 275 :file:`trcsms_my_trc.F90` 276 The routine performs the call to Boundary Conditions and its main purpose is to 277 contain the Source-Minus-Sinks terms due to the biogeochemical processes of the external model. 278 Be aware that lateral boundary conditions are applied in trcnxt routine. 279 280 .. warning:: 281 The routines to compute the light penetration along the water column and 282 the tracer vertical sinking should be defined/called in here, 283 as generalized modules are still missing in the code. 284 285 :file:`trcice_my_trc.F90` 286 Here it is possible to prescribe the tracers concentrations in the sea-ice that 287 will be used as boundary conditions when ice melting occurs (``nn_ice_tr = 1`` in ``&namtrc_ice``). 288 See e.g. the correspondent PISCES subroutine. 289 290 :file:`trcwri_my_trc.F90` 291 This routine performs the output of the model tracers (only those defined in ``&namtrc``) using 292 IOM module (see chapter “Output and Diagnostics” in the :doc:`Reference Manual <cite>`). 293 It is possible to place here the output of additional variables produced by the model, 294 if not done elsewhere in the code, using the call to ``iom_put``. 269 295 270 296 Coupling an external BGC model using NEMO framework … … 273 299 The coupling with an external BGC model through the NEMO compilation framework can be achieved in 274 300 different ways according to the degree of coding complexity of the Biogeochemical model, like e.g., 275 the whole code is made only by one file or it has multiple modules and interfaces spread across several subfolders. 276 277 Beside the 6 core files of MY_TRC module, let’s assume an external BGC model named *MYBGC* and constituted by 278 a rather essential coding structure, likely few Fortran files. 301 the whole code is made only by one file or 302 it has multiple modules and interfaces spread across several subfolders. 303 304 Beside the 6 core files of MY_TRC module, let’s assume an external BGC model named *MYBGC* and 305 constituted by a rather essential coding structure, likely few Fortran files. 279 306 The new coupled configuration name is *NEMO_MYBGC*. 280 307 281 The best solution is to have all files (the modified ``MY_TRC`` routines and the BGC model ones) placed in 282 a unique folder with root ``MYBGCPATH`` and to use the makenemo external readdressing of ``MY_SRC`` folder. 283 284 The coupled configuration listed in ``work_cfgs.txt`` will look like 308 The best solution is to have all files (the modified ``MY_TRC`` routines and the BGC model ones) 309 placed in a unique folder with root ``MYBGCPATH`` and 310 to use the makenemo external readdressing of ``MY_SRC`` folder. 311 312 The coupled configuration listed in :file:`work_cfgs.txt` will look like 285 313 286 314 :: 287 315 288 NEMO_MYBGC OPA_SRC TOP_SRC 316 NEMO_MYBGC OCE TOP 289 317 290 318 and the related ``cpp_MYBGC.fcm`` content will be … … 292 320 .. code-block:: perl 293 321 294 bld::tool::fppkeyskey_iomput key_mpp_mpi key_top295 296 the compilation with ``makenemo`` will be executed through the following syntax322 bld::tool::fppkeys key_iomput key_mpp_mpi key_top 323 324 the compilation with :file:`makenemo` will be executed through the following syntax 297 325 298 326 .. code-block:: console 299 327 300 301 302 The makenemo feature “-e” was introduced to readdress at compilation time the standard MY_SRC folder303 (usually found in NEMO configurations) with a user defined external one. 304 305 The compilation of more articulated BGC model code & infrastructure, like in the case of BFM 306 ([http://www.bfm-community.eu/publications/bfmnemomanual_r1.0_201508.pdf BFM-NEMO coupling manual]),307 requires some additional features.328 $ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>' 329 330 The makenemo feature ``-e`` was introduced to 331 readdress at compilation time the standard MY_SRC folder (usually found in NEMO configurations) with 332 a user defined external one. 333 334 The compilation of more articulated BGC model code & infrastructure, 335 like in the case of BFM (|BFM man|_), requires some additional features. 308 336 309 337 As before, let’s assume a coupled configuration name *NEMO_MYBGC*, 310 but in this case MYBGC model root becomes ``<MYBGCPATH>`` that contains 4 different subfolders for 311 biogeochemistry, named ``initialization``, ``pelagic``, and ``benthic``, and 312 a separate one named ``nemo_coupling`` including the modified ``MY_SRC`` routines. 338 but in this case MYBGC model root becomes :file:`MYBGC` path that 339 contains 4 different subfolders for biogeochemistry, 340 named :file:`initialization`, :file:`pelagic`, and :file:`benthic`, 341 and a separate one named :file:`nemo_coupling` including the modified `MY_SRC` routines. 313 342 The latter folder containing the modified NEMO coupling interface will be still linked using 314 the makenemo “-e”option.343 the makenemo ``-e`` option. 315 344 316 345 In order to include the BGC model subfolders in the compilation of NEMO code, 317 it will be necessary to extend the configuration ``cpp_NEMO_MYBGC.fcm`` file to include the specific paths of 318 ``MYBGC`` folders, as in the following example 346 it will be necessary to extend the configuration :file:`cpp_NEMO_MYBGC.fcm` file to include the specific paths of :file:`MYBGC` folders, as in the following example 319 347 320 348 .. code-block:: perl 321 349 322 323 324 325 326 327 328 329 330 350 bld::tool::fppkeys key_iomput key_mpp_mpi key_top 351 352 src::MYBGC::initialization <MYBGCPATH>/initialization 353 src::MYBGC::pelagic <MYBGCPATH>/pelagic 354 src::MYBGC::benthic <MYBGCPATH>/benthic 355 356 bld::pp::MYBGC 1 357 bld::tool::fppflags::MYBGC %FPPFLAGS 358 bld::tool::fppkeys %bld::tool::fppkeys MYBGC_MACROS 331 359 332 360 where *MYBGC_MACROS* is the space delimited list of macros used in *MYBGC* model for 333 361 selecting/excluding specific parts of the code. 334 The BGC model code will be preprocessed in the configuration ``BLD`` folder as for NEMO,335 but with an independent path, like ``NEMO_MYBGC/BLD/MYBGC/<subforlders>``.362 The BGC model code will be preprocessed in the configuration :file:`BLD` folder as for NEMO, 363 but with an independent path, like :file:`NEMO_MYBGC/BLD/MYBGC/<subforlders>`. 336 364 337 365 The compilation will be performed similarly to in the previous case with the following … … 339 367 .. code-block:: console 340 368 341 $ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>/nemo_coupling' 342 343 Note that, the additional lines specific for the BGC model source and build paths can be written into 344 a separate file, e.g. named ``MYBGC.fcm``, and then simply included in the ``cpp_NEMO_MYBGC.fcm`` as follow 345 346 .. code-block:: perl 347 348 bld::tool::fppkeys key_zdftke key_dynspg_ts key_iomput key_mpp_mpi key_top 349 inc <MYBGCPATH>/MYBGC.fcm 350 351 This will enable a more portable compilation structure for all MYBGC related configurations. 352 353 **Important**: the coupling interface contained in nemo_coupling cannot be added using the FCM syntax, 354 as the same files already exists in NEMO and they are overridden only with the readdressing of MY_SRC contents to 355 avoid compilation conflicts due to duplicate routines. 356 357 All modifications illustrated above, can be easily implemented using shell or python scripting to 358 edit the NEMO configuration CPP.fcm file and to create the BGC model specific FCM compilation file with code paths. 369 $ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>/nemo_coupling' 370 371 .. note:: 372 The additional lines specific for the BGC model source and build paths can be written into 373 a separate file, e.g. named :file:`MYBGC.fcm`, 374 and then simply included in the :file:`cpp_NEMO_MYBGC.fcm` as follow 375 376 .. code-block:: perl 377 378 bld::tool::fppkeys key_zdftke key_dynspg_ts key_iomput key_mpp_mpi key_top 379 inc <MYBGCPATH>/MYBGC.fcm 380 381 This will enable a more portable compilation structure for all MYBGC related configurations. 382 383 .. warning:: 384 The coupling interface contained in :file:`nemo_coupling` cannot be added using the FCM syntax, 385 as the same files already exists in NEMO and they are overridden only with 386 the readdressing of MY_SRC contents to avoid compilation conflicts due to duplicate routines. 387 388 All modifications illustrated above, can be easily implemented using shell or python scripting 389 to edit the NEMO configuration :file:`CPP.fcm` file and 390 to create the BGC model specific FCM compilation file with code paths. 391 392 .. |BFM man| replace:: BFM-NEMO coupling manual -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/TRP/trcnxt.F90
r10425 r12165 139 139 ENDIF 140 140 ! ! Leap-Frog + Asselin filter time stepping 141 IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN ! Euler time-stepping (only swap) 141 IF( (neuler == 0 .AND. kt == nittrc000) ) THEN 142 ! set up for leapfrog on second timestep 143 DO jn = 1, jptra 144 DO jk = 1, jpkm1 145 trb(:,:,jk,jn) = trn(:,:,jk,jn) 146 trn(:,:,jk,jn) = tra(:,:,jk,jn) 147 END DO 148 END DO 149 ELSE IF( ln_top_euler ) THEN 150 ! always doing euler timestepping 142 151 DO jn = 1, jptra 143 152 DO jk = 1, jpkm1 … … 146 155 END DO 147 156 END DO 157 ENDIF 158 IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN ! Euler time-stepping (only swap) 148 159 IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN ! Zero Asselin filter contribution must be explicitly written out since for vvl 149 160 ! ! Asselin filter is output by tra_nxt_vvl that is not called on this time step -
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/trcbdy.F90
r11224 r12165 95 95 END DO 96 96 IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN ! if need to send/recv in at least one direction 97 CALL lbc_lnk( ' bdytra', tsa, 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )97 CALL lbc_lnk( 'trcbdy', tra, 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 98 98 END IF 99 99 !
Note: See TracChangeset
for help on using the changeset viewer.