- Timestamp:
- 2020-09-29T12:41:06+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/r12377_ticket2386
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_adv_pra.F90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r12377_ticket2386
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@13507 sette
-
- Property svn:externals
-
NEMO/branches/2020/r12377_ticket2386/src/ICE/icedyn_adv_pra.F90
r12511 r13540 44 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction 45 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvl , syvl , sxxvl , syyvl , sxyvl ! melt pond lid volume 46 47 47 48 !! * Substitutions … … 55 56 56 57 SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 57 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )58 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 58 59 !!---------------------------------------------------------------------- 59 60 !! ** routine ice_dyn_adv_pra ** … … 81 82 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 82 83 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume 84 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid thickness 83 85 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 84 86 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 85 87 ! 86 INTEGER :: ji, jj, jk, jl, jt! dummy loop indices88 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 87 89 INTEGER :: icycle ! number of sub-timestep for the advection 88 90 REAL(wp) :: zdt ! - - … … 90 92 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 91 93 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 92 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 94 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max 95 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max 96 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: ze_s, zes_max 93 97 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 94 98 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi 95 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp 99 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp, z0vl 96 100 REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es 97 101 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei … … 100 104 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 101 105 ! 102 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 106 ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 107 ! thickness and salinity 108 WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 109 ELSEWHERE ; zs_i(:,:,:) = 0._wp 110 END WHERE 103 111 DO jl = 1, jpl 104 DO_2D _00_00112 DO_2D( 0, 0, 0, 0 ) 105 113 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 106 114 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & … … 115 123 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 116 124 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 125 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 126 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 127 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 128 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 117 129 END_2D 118 130 END DO 119 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 131 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 132 ! 133 ! enthalpies 134 DO jk = 1, nlay_i 135 WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 136 ELSEWHERE ; ze_i(:,:,jk,:) = 0._wp 137 END WHERE 138 END DO 139 DO jk = 1, nlay_s 140 WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 141 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 142 END WHERE 143 END DO 144 DO jl = 1, jpl 145 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 146 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 147 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 148 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 149 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 150 END_3D 151 END DO 152 DO jl = 1, jpl 153 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 154 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 155 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 156 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 157 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 158 END_3D 159 END DO 160 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 161 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 162 ! 120 163 ! 121 164 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! … … 156 199 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 157 200 END DO 158 IF ( ln_pnd_H12 ) THEN 159 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 160 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 201 IF ( ln_pnd_LEV ) THEN 202 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 203 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume 204 IF ( ln_pnd_lids ) THEN 205 z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:) ! Melt pond lid volume 206 ENDIF 161 207 ENDIF 162 208 END DO … … 189 235 END DO 190 236 ! 191 IF ( ln_pnd_ H12) THEN237 IF ( ln_pnd_LEV ) THEN 192 238 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 193 239 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 194 240 CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 195 241 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 242 IF ( ln_pnd_lids ) THEN 243 CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 244 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 245 ENDIF 196 246 ENDIF 197 247 ! !--------------------------------------------! … … 220 270 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 221 271 END DO 222 IF ( ln_pnd_ H12) THEN272 IF ( ln_pnd_LEV ) THEN 223 273 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 224 274 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 225 275 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume 226 276 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 227 ENDIF 277 IF ( ln_pnd_lids ) THEN 278 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 279 CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 280 ENDIF 281 ENDIF 228 282 ! 229 283 ENDIF … … 242 296 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 243 297 END DO 244 IF ( ln_pnd_ H12) THEN298 IF ( ln_pnd_LEV ) THEN 245 299 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 246 300 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 301 IF ( ln_pnd_lids ) THEN 302 pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 303 ENDIF 247 304 ENDIF 248 305 END DO … … 250 307 ! derive open water from ice concentration 251 308 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 252 DO_2D _00_00309 DO_2D( 0, 0, 0, 0 ) 253 310 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 254 311 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 255 312 END_2D 256 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. )313 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp ) 257 314 ! 258 315 ! --- Ensure non-negative fields --- ! 259 316 ! Remove negative values (conservation is ensured) 260 317 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 261 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, p e_s, pe_i )318 CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 262 319 ! 263 320 ! --- Make sure ice thickness is not too big --- ! 264 321 ! (because ice thickness can be too large where ice concentration is very small) 265 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 322 CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 323 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 266 324 ! 267 325 ! --- Ensure snow load is not too big --- ! … … 305 363 ! 306 364 ! Limitation of moments. 307 DO_2D _00_11365 DO_2D( 0, 0, 1, 1 ) 308 366 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 309 367 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) … … 325 383 326 384 ! Calculate fluxes and moments between boxes i<-->i+1 327 DO_2D _00_11385 DO_2D( 0, 0, 1, 1 ) ! Flux from i to i+1 WHEN u GT 0 328 386 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 329 387 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) … … 350 408 END_2D 351 409 352 DO_2D _00_10410 DO_2D( 0, 0, 1, 0 ) ! Flux from i+1 to i when u LT 0. 353 411 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 354 412 zalg (ji,jj) = zalf … … 369 427 END_2D 370 428 371 DO_2D _00_00429 DO_2D( 0, 0, 0, 0 ) ! Readjust moments remaining in the box. 372 430 zbt = zbet(ji-1,jj) 373 431 zbt1 = 1.0 - zbet(ji-1,jj) … … 383 441 384 442 ! Put the temporary moments into appropriate neighboring boxes. 385 DO_2D _00_00443 DO_2D( 0, 0, 0, 0 ) ! Flux from i to i+1 IF u GT 0. 386 444 zbt = zbet(ji-1,jj) 387 445 zbt1 = 1.0 - zbet(ji-1,jj) … … 403 461 END_2D 404 462 405 DO_2D _00_00463 DO_2D( 0, 0, 0, 0 ) ! Flux from i+1 to i IF u LT 0. 406 464 zbt = zbet(ji,jj) 407 465 zbt1 = 1.0 - zbet(ji,jj) … … 425 483 426 484 !-- Lateral boundary conditions 427 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1. , ps0 , 'T', 1.&428 & , psx , 'T', -1. , psy , 'T', -1.& ! caution gradient ==> the sign changes429 & , psxx , 'T', 1. , psyy, 'T', 1. , psxy, 'T', 1.)485 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 486 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 487 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 430 488 ! 431 489 END SUBROUTINE adv_x … … 462 520 ! 463 521 ! Limitation of moments. 464 DO_2D _11_00522 DO_2D( 1, 1, 0, 0 ) 465 523 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 466 524 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) … … 482 540 483 541 ! Calculate fluxes and moments between boxes j<-->j+1 484 DO_2D _11_00542 DO_2D( 1, 1, 0, 0 ) ! Flux from j to j+1 WHEN v GT 0 485 543 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 486 544 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) … … 507 565 END_2D 508 566 ! 509 DO_2D _10_00567 DO_2D( 1, 0, 0, 0 ) ! Flux from j+1 to j when v LT 0. 510 568 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 511 569 zalg (ji,jj) = zalf … … 527 585 528 586 ! Readjust moments remaining in the box. 529 DO_2D _00_00587 DO_2D( 0, 0, 0, 0 ) 530 588 zbt = zbet(ji,jj-1) 531 589 zbt1 = ( 1.0 - zbet(ji,jj-1) ) … … 541 599 542 600 ! Put the temporary moments into appropriate neighboring boxes. 543 DO_2D _00_00601 DO_2D( 0, 0, 0, 0 ) ! Flux from j to j+1 IF v GT 0. 544 602 zbt = zbet(ji,jj-1) 545 603 zbt1 = 1.0 - zbet(ji,jj-1) … … 562 620 END_2D 563 621 564 DO_2D _00_00622 DO_2D( 0, 0, 0, 0 ) ! Flux from j+1 to j IF v LT 0. 565 623 zbt = zbet(ji,jj) 566 624 zbt1 = 1.0 - zbet(ji,jj) … … 584 642 585 643 !-- Lateral boundary conditions 586 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1. , ps0 , 'T', 1.&587 & , psx , 'T', -1. , psy , 'T', -1.& ! caution gradient ==> the sign changes588 & , psxx , 'T', 1. , psyy, 'T', 1. , psxy, 'T', 1.)644 CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1.0_wp, ps0 , 'T', 1.0_wp & 645 & , psx , 'T', -1.0_wp, psy , 'T', -1.0_wp & ! caution gradient ==> the sign changes 646 & , psxx , 'T', 1.0_wp, psyy, 'T', 1.0_wp , psxy, 'T', 1.0_wp ) 589 647 ! 590 648 END SUBROUTINE adv_y 591 649 592 650 593 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 651 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 652 & pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 594 653 !!------------------------------------------------------------------- 595 654 !! *** ROUTINE Hbig *** … … 605 664 !! ** input : Max thickness of the surrounding 9-points 606 665 !!------------------------------------------------------------------- 607 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 608 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 609 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip 666 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 667 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts 668 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max 669 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max 670 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 610 671 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 611 ! 612 INTEGER :: ji, jj, jl ! dummy loop indices 613 REAL(wp) :: z1_dt, zhip, zhi, zhs, zfra 672 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 673 ! 674 INTEGER :: ji, jj, jk, jl ! dummy loop indices 675 REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 614 676 !!------------------------------------------------------------------- 615 677 ! … … 617 679 ! 618 680 DO jl = 1, jpl 619 620 DO_2D_11_11 681 DO_2D( 1, 1, 1, 1 ) 621 682 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 622 683 ! 623 684 ! ! -- check h_ip -- ! 624 685 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 625 IF( ln_pnd_ H12.AND. pv_ip(ji,jj,jl) > 0._wp ) THEN686 IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 626 687 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 627 688 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN … … 650 711 ENDIF 651 712 ! 713 ! ! -- check s_i -- ! 714 ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 715 zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 716 IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 717 zfra = psi_max(ji,jj,jl) / zsi 718 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 719 psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 720 ENDIF 721 ! 652 722 ENDIF 653 723 END_2D 654 724 END DO 725 ! 726 ! ! -- check e_i/v_i -- ! 727 DO jl = 1, jpl 728 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 729 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 730 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 731 zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 732 IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 733 zfra = pei_max(ji,jj,jk,jl) / zei 734 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 735 pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 736 ENDIF 737 ENDIF 738 END_3D 739 END DO 740 ! ! -- check e_s/v_s -- ! 741 DO jl = 1, jpl 742 DO_3D( 1, 1, 1, 1, 1, nlay_s ) 743 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 744 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 745 zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 746 IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 747 zfra = pes_max(ji,jj,jk,jl) / zes 748 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 749 pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 750 ENDIF 751 ENDIF 752 END_3D 753 END DO 655 754 ! 656 755 END SUBROUTINE Hbig … … 684 783 ! -- check snow load -- ! 685 784 DO jl = 1, jpl 686 DO_2D _11_11785 DO_2D( 1, 1, 1, 1 ) 687 786 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 688 787 ! … … 724 823 & sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) , & 725 824 & sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) , & 726 & sxap(jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , & 727 & sxvp(jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , & 825 & sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , & 826 & sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , & 827 & sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) , & 728 828 ! 729 829 & sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & … … 772 872 ! 773 873 ! ! ice thickness 774 CALL iom_get( numrir, jpdom_auto glo, 'sxice' , sxice )775 CALL iom_get( numrir, jpdom_auto glo, 'syice' , syice )776 CALL iom_get( numrir, jpdom_auto glo, 'sxxice', sxxice )777 CALL iom_get( numrir, jpdom_auto glo, 'syyice', syyice )778 CALL iom_get( numrir, jpdom_auto glo, 'sxyice', sxyice )874 CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice ) 875 CALL iom_get( numrir, jpdom_auto, 'syice' , syice ) 876 CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice ) 877 CALL iom_get( numrir, jpdom_auto, 'syyice', syyice ) 878 CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice ) 779 879 ! ! snow thickness 780 CALL iom_get( numrir, jpdom_auto glo, 'sxsn' , sxsn )781 CALL iom_get( numrir, jpdom_auto glo, 'sysn' , sysn )782 CALL iom_get( numrir, jpdom_auto glo, 'sxxsn' , sxxsn )783 CALL iom_get( numrir, jpdom_auto glo, 'syysn' , syysn )784 CALL iom_get( numrir, jpdom_auto glo, 'sxysn' , sxysn )880 CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn ) 881 CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn ) 882 CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn ) 883 CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn ) 884 CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn ) 785 885 ! ! ice concentration 786 CALL iom_get( numrir, jpdom_auto glo, 'sxa' , sxa )787 CALL iom_get( numrir, jpdom_auto glo, 'sya' , sya )788 CALL iom_get( numrir, jpdom_auto glo, 'sxxa' , sxxa )789 CALL iom_get( numrir, jpdom_auto glo, 'syya' , syya )790 CALL iom_get( numrir, jpdom_auto glo, 'sxya' , sxya )886 CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa ) 887 CALL iom_get( numrir, jpdom_auto, 'sya' , sya ) 888 CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa ) 889 CALL iom_get( numrir, jpdom_auto, 'syya' , syya ) 890 CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya ) 791 891 ! ! ice salinity 792 CALL iom_get( numrir, jpdom_auto glo, 'sxsal' , sxsal )793 CALL iom_get( numrir, jpdom_auto glo, 'sysal' , sysal )794 CALL iom_get( numrir, jpdom_auto glo, 'sxxsal', sxxsal )795 CALL iom_get( numrir, jpdom_auto glo, 'syysal', syysal )796 CALL iom_get( numrir, jpdom_auto glo, 'sxysal', sxysal )892 CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal ) 893 CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal ) 894 CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal ) 895 CALL iom_get( numrir, jpdom_auto, 'syysal', syysal ) 896 CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal ) 797 897 ! ! ice age 798 CALL iom_get( numrir, jpdom_auto glo, 'sxage' , sxage )799 CALL iom_get( numrir, jpdom_auto glo, 'syage' , syage )800 CALL iom_get( numrir, jpdom_auto glo, 'sxxage', sxxage )801 CALL iom_get( numrir, jpdom_auto glo, 'syyage', syyage )802 CALL iom_get( numrir, jpdom_auto glo, 'sxyage', sxyage )898 CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage ) 899 CALL iom_get( numrir, jpdom_auto, 'syage' , syage ) 900 CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage ) 901 CALL iom_get( numrir, jpdom_auto, 'syyage', syyage ) 902 CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage ) 803 903 ! ! snow layers heat content 804 904 DO jk = 1, nlay_s 805 905 WRITE(zchar1,'(I2.2)') jk 806 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)807 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)808 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:)809 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:)810 znam = 'sxyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:)906 znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:) 907 znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:) 908 znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:) 909 znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:) 910 znam = 'sxyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:) 811 911 END DO 812 912 ! ! ice layers heat content 813 913 DO jk = 1, nlay_i 814 914 WRITE(zchar1,'(I2.2)') jk 815 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)816 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)817 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:)818 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:)819 znam = 'sxye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto glo, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:)915 znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:) 916 znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:) 917 znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:) 918 znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:) 919 znam = 'sxye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:) 820 920 END DO 821 921 ! 822 IF( ln_pnd_H12 ) THEN ! melt pond fraction 823 CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap ) 824 CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap ) 825 CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap ) 826 CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap ) 827 CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap ) 828 ! ! melt pond volume 829 CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp ) 830 CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp ) 831 CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp ) 832 CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 833 CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 922 IF( ln_pnd_LEV ) THEN ! melt pond fraction 923 IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 924 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap ) 925 CALL iom_get( numrir, jpdom_auto, 'syap' , syap ) 926 CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 927 CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 928 CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 929 ! ! melt pond volume 930 CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp ) 931 CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp ) 932 CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 933 CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 934 CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 935 ELSE 936 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 937 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 938 ENDIF 939 ! 940 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 941 IF( iom_varid( numror, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 942 CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl ) 943 CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl ) 944 CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 945 CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) 946 CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl ) 947 ELSE 948 sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume 949 ENDIF 950 ENDIF 834 951 ENDIF 835 952 ! … … 845 962 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 846 963 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content 847 IF( ln_pnd_H12 ) THEN 848 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 849 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 964 IF( ln_pnd_LEV ) THEN 965 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 966 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume 967 IF ( ln_pnd_lids ) THEN 968 sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume 969 ENDIF 850 970 ENDIF 851 971 ENDIF … … 910 1030 END DO 911 1031 ! 912 IF( ln_pnd_ H12) THEN ! melt pond fraction1032 IF( ln_pnd_LEV ) THEN ! melt pond fraction 913 1033 CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap ) 914 1034 CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap ) … … 922 1042 CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 923 1043 CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 1044 ! 1045 IF ( ln_pnd_lids ) THEN ! melt pond lid volume 1046 CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl ) 1047 CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl ) 1048 CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl ) 1049 CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl ) 1050 CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl ) 1051 ENDIF 924 1052 ENDIF 925 1053 !
Note: See TracChangeset
for help on using the changeset viewer.
