- Timestamp:
- 2019-11-26T12:08:01+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
- Files:
-
- 20 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cpl_oasis3.F90
r10582 r11963 114 114 !------------------------------------------------------------------ 115 115 CALL oasis_init_comp ( ncomp_id, TRIM(cd_modname), nerror ) 116 IF 116 IF( nerror /= OASIS_Ok ) & 117 117 CALL oasis_abort (ncomp_id, 'cpl_init', 'Failure in oasis_init_comp') 118 118 … … 122 122 123 123 CALL oasis_get_localcomm ( kl_comm, nerror ) 124 IF 124 IF( nerror /= OASIS_Ok ) & 125 125 CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' ) 126 126 ! … … 149 149 150 150 ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 151 IF 151 IF( ltmp_wapatch ) THEN 152 152 nldi_save = nldi ; nlei_save = nlei 153 153 nldj_save = nldj ; nlej_save = nlej … … 217 217 ! 218 218 DO ji = 1, ksnd 219 IF 219 IF( ssnd(ji)%laction ) THEN 220 220 221 221 IF( ssnd(ji)%nct > nmaxcat ) THEN … … 228 228 DO jm = 1, kcplmodel 229 229 230 IF 230 IF( ssnd(ji)%nct .GT. 1 ) THEN 231 231 WRITE(cli2,'(i2.2)') jc 232 232 zclname = TRIM(ssnd(ji)%clname)//'_cat'//cli2 … … 234 234 zclname = ssnd(ji)%clname 235 235 ENDIF 236 IF 236 IF( kcplmodel > 1 ) THEN 237 237 WRITE(cli2,'(i2.2)') jm 238 238 zclname = 'model'//cli2//'_'//TRIM(zclname) … … 241 241 IF( agrif_fixed() /= 0 ) THEN 242 242 zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname) 243 END 243 ENDIF 244 244 #endif 245 245 IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out 246 246 CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), & 247 247 & OASIS_Out , ishape , OASIS_REAL, nerror ) 248 IF 248 IF( nerror /= OASIS_Ok ) THEN 249 249 WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname) 250 250 CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' ) … … 262 262 ! 263 263 DO ji = 1, krcv 264 IF 264 IF( srcv(ji)%laction ) THEN 265 265 266 266 IF( srcv(ji)%nct > nmaxcat ) THEN … … 273 273 DO jm = 1, kcplmodel 274 274 275 IF 275 IF( srcv(ji)%nct .GT. 1 ) THEN 276 276 WRITE(cli2,'(i2.2)') jc 277 277 zclname = TRIM(srcv(ji)%clname)//'_cat'//cli2 … … 279 279 zclname = srcv(ji)%clname 280 280 ENDIF 281 IF 281 IF( kcplmodel > 1 ) THEN 282 282 WRITE(cli2,'(i2.2)') jm 283 283 zclname = 'model'//cli2//'_'//TRIM(zclname) … … 286 286 IF( agrif_fixed() /= 0 ) THEN 287 287 zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname) 288 END 288 ENDIF 289 289 #endif 290 290 IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In 291 291 CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part , (/ 2, 1 /), & 292 292 & OASIS_In , ishape , OASIS_REAL, nerror ) 293 IF 293 IF( nerror /= OASIS_Ok ) THEN 294 294 WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname) 295 295 CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' ) … … 310 310 IF( nerror /= OASIS_Ok ) CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef') 311 311 ! 312 IF 312 IF( ltmp_wapatch ) THEN 313 313 nldi = nldi_save ; nlei = nlei_save 314 314 nldj = nldj_save ; nlej = nlej_save … … 332 332 !!-------------------------------------------------------------------- 333 333 ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 334 IF 334 IF( ltmp_wapatch ) THEN 335 335 nldi_save = nldi ; nlei_save = nlei 336 336 nldj_save = nldj ; nlej_save = nlej … … 349 349 CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo ) 350 350 351 IF 352 IF 351 IF( ln_ctl ) THEN 352 IF( kinfo == OASIS_Sent .OR. kinfo == OASIS_ToRest .OR. & 353 353 & kinfo == OASIS_SentOut .OR. kinfo == OASIS_ToRestOut ) THEN 354 354 WRITE(numout,*) '****************' … … 368 368 ENDDO 369 369 ENDDO 370 IF 370 IF( ltmp_wapatch ) THEN 371 371 nldi = nldi_save ; nlei = nlei_save 372 372 nldj = nldj_save ; nlej = nlej_save … … 393 393 !!-------------------------------------------------------------------- 394 394 ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 395 IF 395 IF( ltmp_wapatch ) THEN 396 396 nldi_save = nldi ; nlei_save = nlei 397 397 nldj_save = nldj ; nlej_save = nlej … … 403 403 ! 404 404 DO jc = 1, srcv(kid)%nct 405 IF 405 IF( ltmp_wapatch ) THEN 406 406 IF( nimpp == 1 ) nldi = 1 407 407 IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi … … 420 420 & kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut 421 421 422 IF 422 IF( ln_ctl ) WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm) 423 423 424 IF 424 IF( llaction ) THEN 425 425 426 426 kinfo = OASIS_Rcv … … 432 432 ENDIF 433 433 434 IF 434 IF( ln_ctl ) THEN 435 435 WRITE(numout,*) '****************' 436 436 WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname … … 450 450 ENDDO 451 451 452 IF 452 IF( ltmp_wapatch ) THEN 453 453 nldi = nldi_save ; nlei = nlei_save 454 454 nldj = nldj_save ; nlej = nlej_save … … 483 483 ! 484 484 DO ji = 1, nsnd 485 IF 485 IF(ssnd(ji)%laction ) THEN 486 486 DO jm = 1, ncplmodel 487 487 IF( ssnd(ji)%nid(1,jm) /= -1 ) THEN … … 495 495 ENDDO 496 496 DO ji = 1, nrcv 497 IF 497 IF(srcv(ji)%laction ) THEN 498 498 DO jm = 1, ncplmodel 499 499 IF( srcv(ji)%nid(1,jm) /= -1 ) THEN … … 529 529 ! 530 530 DEALLOCATE( exfld ) 531 IF 531 IF(nstop == 0) THEN 532 532 CALL oasis_terminate( nerror ) 533 533 ELSE -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cyclone.F90
r10068 r11963 137 137 zhemi = SIGN( 1. , zrlat ) 138 138 zinfl = 15.* rad ! clim inflow angle in Tropical Cyclones 139 IF 139 IF( vortex == 0 ) THEN 140 140 141 141 ! Vortex Holland reconstruct wind at each lon-lat position … … 157 157 & + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 158 158 159 IF 159 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 160 160 ! shape of the wind profile 161 161 zztmp = ( zrmw / ( zdist + 1.e-12 ) )**zb 162 162 zztmp = zvmax * SQRT( zztmp * EXP(1. - zztmp) ) 163 163 164 IF 164 IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 165 165 zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 166 166 ENDIF 167 167 168 168 ! !!! KILL EQ WINDS 169 ! IF 169 ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 170 170 ! zztmp = 0. ! winds in other hemisphere 171 ! IF 172 ! ENDIF 173 ! IF 171 ! IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0. ! kill between 5N-5S 172 ! ENDIF 173 ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 174 174 ! zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) ) 175 175 ! !linear to zero between 10 and 5 … … 177 177 ! !!! / KILL EQ 178 178 179 IF 179 IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 180 180 181 181 zwnd_t = COS( zinfl ) * zztmp … … 196 196 END DO 197 197 198 ELSE IF 198 ELSE IF( vortex == 1 ) THEN 199 199 200 200 ! Vortex Willoughby reconstruct wind at each lon-lat position … … 206 206 zn = 2.1340 + 0.0077*zvmax - 0.4522*LOG(zrmw/1000.) - 0.0038*ABS( ztct(jtc,jp_lat) ) 207 207 zA = 0.5913 + 0.0029*zvmax - 0.1361*LOG(zrmw/1000.) - 0.0042*ABS( ztct(jtc,jp_lat) ) 208 IF 208 IF(zA < 0) THEN 209 209 zA=0 210 210 ENDIF … … 218 218 & + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 219 219 220 IF 220 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 221 221 222 222 ! shape of the wind profile 223 IF 223 IF(zdist <= zrmw) THEN ! inside the Radius of Maximum Wind 224 224 zztmp = zvmax * (zdist/zrmw)**zn 225 225 ELSE … … 227 227 ENDIF 228 228 229 IF 229 IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 230 230 zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 231 231 ENDIF 232 232 233 233 ! !!! KILL EQ WINDS 234 ! IF 234 ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 235 235 ! zztmp = 0. ! winds in other hemisphere 236 ! IF 237 ! ENDIF 238 ! IF 236 ! IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0. ! kill between 5N-5S 237 ! ENDIF 238 ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 239 239 ! zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) ) 240 240 ! !linear to zero between 10 and 5 … … 242 242 ! !!! / KILL EQ 243 243 244 IF 244 IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 245 245 246 246 zwnd_t = COS( zinfl ) * zztmp -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/fldread.F90
r11831 r11963 167 167 IF( PRESENT(kit) ) ll_firstcall = ll_firstcall .and. kit == 1 168 168 169 IF 169 IF( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc 170 170 ELSE ; it_offset = 0 171 171 ENDIF … … 389 389 ENDIF 390 390 ! 391 IF 391 IF( sdjf%cltype(1:4) == 'week' ) THEN 392 392 isec_week = isec_week + ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week 393 393 llprevmth = isec_week > nsec_month ! longer time since the beginning of the week than the month … … 464 464 ENDIF 465 465 ! 466 IF 466 IF( nn_components == jp_iam_sas ) THEN ; it_offset = nn_fsbc 467 467 ELSE ; it_offset = 0 468 468 ENDIF … … 656 656 ENDIF 657 657 CASE DEFAULT 658 IF 658 IF(lk_c1d .AND. lmoor ) THEN 659 659 IF( sdjf%ln_tint ) THEN 660 660 CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) ) … … 1071 1071 imonth = kmonth 1072 1072 iday = kday 1073 IF 1073 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1074 1074 isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 ) 1075 1075 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month … … 1080 1080 ENDIF 1081 1081 ELSE ! use current day values 1082 IF 1082 IF( sdjf%cltype(1:4) == 'week' ) THEN ! find the day of the beginning of the week 1083 1083 isec_week = ksec_week( sdjf%cltype(6:8) ) ! second since the beginning of the week 1084 1084 llprevmth = isec_week > nsec_month ! longer time since beginning of the week than the month … … 1318 1318 1319 1319 !! get dimensions 1320 IF 1320 IF( SIZE(sd%fnow, 3) > 1 ) THEN 1321 1321 ALLOCATE( ddims(4) ) 1322 1322 ELSE … … 1331 1331 1332 1332 CALL iom_open ( sd%wgtname, inum ) ! interpolation weights 1333 IF 1333 IF( inum > 0 ) THEN 1334 1334 1335 1335 !! determine whether we have an east-west cyclic grid … … 1663 1663 END DO 1664 1664 1665 IF 1665 IF(ref_wgts(kw)%numwgt .EQ. 16) THEN 1666 1666 1667 1667 !! fix up halo points that we couldnt read from file … … 1746 1746 END DO 1747 1747 ! 1748 END 1748 ENDIF 1749 1749 ! 1750 1750 END SUBROUTINE fld_interp -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcapr.F90
r11831 r11963 103 103 ! 104 104 ! !* control check 105 IF 105 IF( ln_apr_obc ) THEN 106 106 IF(lwp) WRITE(numout,*) ' Inverse barometer added to OBC ssh data' 107 107 ENDIF -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
r11962 r11963 259 259 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 260 260 ! 261 IF 261 IF( ln_wave ) THEN 262 262 !Activated wave module but neither drag nor stokes drift activated 263 IF 263 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 264 264 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 265 265 !drag coefficient read from wave model definable only with mfs bulk formulae and core 266 ELSEIF 266 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 267 267 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 268 ELSEIF 268 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 269 269 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 270 270 ENDIF 271 271 ELSE 272 IF 272 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 273 273 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 274 274 & 'with drag coefficient (ln_cdgw =T) ' , & … … 481 481 zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 482 482 483 IF 483 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 484 484 zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action 485 485 zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action 486 END 486 ENDIF 487 487 488 488 !LB: … … 492 492 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 493 493 CASE( np_humi_dpt ) 494 !IF 494 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 495 495 zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 496 496 CASE( np_humi_rlh ) 497 !IF 497 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 498 498 zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file 499 499 END SELECT … … 506 506 507 507 508 IF 508 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 509 509 510 510 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point … … 543 543 tsk(:,:) = zst(:,:) 544 544 545 ELSE !IF 545 ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 546 546 547 547 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point … … 567 567 END SELECT 568 568 569 END IF ! IF( ln_skin_cs .OR. ln_skin_wl )569 ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 570 570 571 571 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. … … 576 576 t_zu(:,:) = ztpot(:,:) 577 577 q_zu(:,:) = zqair(:,:) 578 END 578 ENDIF 579 579 580 580 … … 660 660 #endif 661 661 ! 662 !!#LB: NO WHY???? IF 662 !!#LB: NO WHY???? IF( nn_ice == 0 ) THEN 663 663 CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB 664 664 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean … … 674 674 CALL iom_put( 'snowpre', sprecip ) ! Snow 675 675 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 676 IF 676 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 677 677 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 678 678 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 679 END 679 ENDIF 680 680 !!#LB. ENDIF 681 681 ! … … 819 819 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 820 820 CASE( np_humi_dpt ) 821 !IF 821 !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm 822 822 zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 823 823 CASE( np_humi_rlh ) 824 !IF 824 !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm 825 825 zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file 826 826 END SELECT -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90
r11962 r11963 68 68 INTEGER :: ierr 69 69 !!--------------------------------------------------------------------- 70 IF 70 IF( l_use_wl ) THEN 71 71 ierr = 0 72 72 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) … … 76 76 dT_wl(:,:) = 0._wp 77 77 Hz_wl(:,:) = Hwl_max 78 END 79 IF 78 ENDIF 79 IF( l_use_cs ) THEN 80 80 ierr = 0 81 81 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 82 82 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P0_INIT => allocation of dT_cs failed!' ) 83 83 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 84 END 84 ENDIF 85 85 END SUBROUTINE sbcblk_algo_coare3p0_init 86 86 … … 190 190 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0' 191 191 !!---------------------------------------------------------------------------------- 192 IF 192 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl) 193 193 194 194 l_zt_equal_zu = .FALSE. … … 197 197 198 198 !! Initializations for cool skin and warm layer: 199 IF 199 IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 200 200 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 201 201 202 IF 202 IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 203 203 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 204 204 205 IF 205 IF( l_use_cs .OR. l_use_wl ) THEN 206 206 ALLOCATE ( zsst(jpi,jpj) ) 207 207 zsst = T_s ! backing up the bulk SST 208 208 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 209 209 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 210 END 210 ENDIF 211 211 212 212 … … 265 265 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 266 266 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 267 END 267 ENDIF 268 268 269 269 !! ITERATION BLOCK … … 288 288 zeta_t = zt*ztmp0 289 289 zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 290 END 290 ENDIF 291 291 292 292 !! Adjustment the wind at 10m (not needed in the current algo form): 293 !IF 293 !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 294 294 295 295 !! Roughness lengthes z0, z0t (z0q = z0t) : … … 299 299 300 300 ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 301 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #L OLO: some use 1.15 not 1.1 !!!301 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LB: some use 1.15 not 1.1 !!! 302 302 z0t = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 303 303 … … 315 315 t_zu = t_zt - t_star/vkarmn*ztmp1 316 316 q_zu = q_zt - q_star/vkarmn*ztmp1 317 END 317 ENDIF 318 318 319 319 … … 330 330 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 331 331 332 END 332 ENDIF 333 333 334 334 IF( l_use_wl ) THEN … … 343 343 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 344 344 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 345 END 345 ENDIF 346 346 347 347 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 348 348 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 349 349 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 350 END 350 ENDIF 351 351 352 352 END DO !DO j_itt = 1, nb_itt … … 365 365 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 366 366 367 IF 368 IF 369 IF 370 371 IF 367 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 368 IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 369 IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 370 371 IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 372 372 373 373 END SUBROUTINE turb_coare3p0 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90
r11962 r11963 68 68 INTEGER :: ierr 69 69 !!--------------------------------------------------------------------- 70 IF 70 IF( l_use_wl ) THEN 71 71 ierr = 0 72 72 ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) … … 76 76 dT_wl(:,:) = 0._wp 77 77 Hz_wl(:,:) = Hwl_max 78 END 79 IF 78 ENDIF 79 IF( l_use_cs ) THEN 80 80 ierr = 0 81 81 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 82 82 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P6_INIT => allocation of dT_cs failed!' ) 83 83 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 84 END 84 ENDIF 85 85 END SUBROUTINE sbcblk_algo_coare3p6_init 86 86 … … 190 190 CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6' 191 191 !!---------------------------------------------------------------------------------- 192 IF 192 IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 193 193 194 194 l_zt_equal_zu = .FALSE. … … 197 197 198 198 !! Initializations for cool skin and warm layer: 199 IF 199 IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 200 200 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 201 201 202 IF 202 IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 203 203 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 204 204 205 IF 205 IF( l_use_cs .OR. l_use_wl ) THEN 206 206 ALLOCATE ( zsst(jpi,jpj) ) 207 207 zsst = T_s ! backing up the bulk SST 208 208 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 209 209 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 210 END 210 ENDIF 211 211 212 212 … … 265 265 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 266 266 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 267 END 267 ENDIF 268 268 269 269 !! ITERATION BLOCK … … 288 288 zeta_t = zt*ztmp0 289 289 zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 290 END 290 ENDIF 291 291 292 292 !! Adjustment the wind at 10m (not needed in the current algo form): 293 !IF 293 !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 294 294 295 295 !! Roughness lengthes z0, z0t (z0q = z0t) : … … 315 315 t_zu = t_zt - t_star/vkarmn*ztmp1 316 316 q_zu = q_zt - q_star/vkarmn*ztmp1 317 END 317 ENDIF 318 318 319 319 … … 330 330 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 331 331 332 END 332 ENDIF 333 333 334 334 IF( l_use_wl ) THEN … … 343 343 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 344 344 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 345 END 345 ENDIF 346 346 347 347 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 348 348 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 349 349 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 350 END 350 ENDIF 351 351 352 352 END DO !DO j_itt = 1, nb_itt … … 365 365 IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 366 366 367 IF 368 IF 369 IF 370 371 IF 367 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 368 IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 369 IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 370 371 IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 372 372 373 373 END SUBROUTINE turb_coare3p6 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r11962 r11963 74 74 INTEGER :: ierr 75 75 !!--------------------------------------------------------------------- 76 IF 76 IF( l_use_wl ) THEN 77 77 ierr = 0 78 78 ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) … … 80 80 dT_wl(:,:) = 0._wp 81 81 Hz_wl(:,:) = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 82 END 83 IF 82 ENDIF 83 IF( l_use_cs ) THEN 84 84 ierr = 0 85 85 ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 86 86 IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 87 87 dT_cs(:,:) = -0.25_wp ! First guess of skin correction 88 END 88 ENDIF 89 89 END SUBROUTINE sbcblk_algo_ecmwf_init 90 90 … … 195 195 !!---------------------------------------------------------------------------------- 196 196 197 IF 197 IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 198 198 199 199 l_zt_equal_zu = .FALSE. … … 201 201 202 202 !! Initializations for cool skin and warm layer: 203 IF 203 IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 204 204 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 205 205 206 IF 206 IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 207 207 & CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 208 208 209 IF 209 IF( l_use_cs .OR. l_use_wl ) THEN 210 210 ALLOCATE ( zsst(jpi,jpj) ) 211 211 zsst = T_s ! backing up the bulk SST 212 212 IF( l_use_cs ) T_s = T_s - 0.25_wp ! First guess of correction 213 213 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 214 END 214 ENDIF 215 215 216 216 … … 270 270 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 271 271 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 272 END 272 ENDIF 273 273 274 274 … … 293 293 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 294 294 !! Note: it is slightly different that the L we would get with the usual 295 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO295 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) 296 296 297 297 !! Update func_m with new Linv: … … 335 335 ztmp1 = LOG(zt/zu) + ztmp2 336 336 q_zu = q_zt - q_star/vkarmn*ztmp1 337 END 337 ENDIF 338 338 339 339 !! Updating because of updated z0 and z0t and new Linv... … … 355 355 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 356 356 357 END 357 ENDIF 358 358 359 359 IF( l_use_wl ) THEN … … 366 366 IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 367 367 q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 368 END 368 ENDIF 369 369 370 370 IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 371 371 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 372 372 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 373 END 373 ENDIF 374 374 375 375 END DO !DO j_itt = 1, nb_itt … … 384 384 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 385 385 386 IF 387 IF 388 IF 389 390 IF 386 IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 387 IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 388 IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 389 390 IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 391 391 392 392 END SUBROUTINE turb_ecmwf -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
r11845 r11963 171 171 q_zu = q_zt - ztmp2/vkarmn*stab ! ztmp2 is still q* L&Y 2004 eq.(9c) 172 172 q_zu = max(0._wp, q_zu) 173 END 173 ENDIF 174 174 175 175 ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
r11962 r11963 541 541 pQns(ji,jj) = zQlat + zQsen + zQlw 542 542 543 IF 543 IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 544 544 END DO 545 545 END DO … … 594 594 pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 595 595 596 IF 597 IF 596 IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 597 IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 598 598 599 599 END DO … … 647 647 pQlat = L_vap(pTs) * zevap 648 648 649 IF 650 IF 649 IF( PRESENT(pEvap) ) pEvap = - zevap 650 IF( PRESENT(prhoa) ) prhoa = zrho 651 651 652 652 END SUBROUTINE BULK_FORMULA_SCLR -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
r11962 r11963 2 2 !!====================================================================== 3 3 !! *** MODULE sbcblk_skin_coare *** 4 !! Computes:5 !! * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer6 !! scheme used at ECMWF7 !! Using formulation/param. of COARE 3.6 (Fairall et al., 2019)8 4 !! 9 !! From routine turb_ecmwf maintained and developed in AeroBulk10 !! (https://github.com/brodeau/aerobulk)5 !! Module that gathers the cool-skin and warm-layer parameterization used 6 !! in the COARE family of bulk parameterizations. 11 7 !! 12 !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 8 !! Based on the last update for version COARE 3.6 (Fairall et al., 2019) 9 !! 10 !! Module 'sbcblk_skin_coare' also maintained and developed in AeroBulk (as 11 !! 'mod_skin_coare') 12 !! (https://github.com/brodeau/aerobulk) !! 13 !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 13 14 !!---------------------------------------------------------------------- 14 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 15 !!---------------------------------------------------------------------- 16 17 !!---------------------------------------------------------------------- 18 !! cswl_ecmwf : computes the surface skin temperature (aka SSST) 19 !! based on the cool-skin/warm-layer scheme used at ECMWF 15 !! History : 4.x ! 2019-11 (L.Brodeau) Original code 20 16 !!---------------------------------------------------------------------- 21 17 USE oce ! ocean dynamics and tracers … … 24 20 USE sbc_oce ! Surface boundary condition: ocean fields 25 21 26 USE sbcblk_phy ! LOLO: all thermodynamics functions, rho_air, q_sat, etc...27 28 USE sbcdcy ! LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE)29 22 USE sbcblk_phy ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc) 23 24 USE sbcdcy !#LB: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE) 25 30 26 USE lib_mpp ! distribued memory computing library 31 27 USE in_out_manager ! I/O manager … … 38 34 39 35 !! Cool-skin related parameters: 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 41 & dT_cs !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect 37 ! ! => temperature difference between air-sea interface (z=0) 38 ! ! and right below viscous layer (z=delta) 42 39 43 40 !! Warm-layer related parameters: 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 45 & dT_wl, & !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 46 & Hz_wl, & !: depth of warm-layer [m] 47 & Qnt_ac, & !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) 48 & Tau_ac !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect 42 ! ! => difference between "almost surface (right below 43 ! ! viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m] 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Qnt_ac !: time integral / accumulated heat stored by the warm layer 46 ! ! Qxdt => [J/m^2] (reset to zero every midnight) 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Tau_ac !: time integral / accumulated momentum 48 ! ! Tauxdt => [N.s/m^2] (reset to zero every midnight) 49 49 50 50 REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp !: maximum depth of warm layer (adjustable) … … 85 85 !!--------------------------------------------------------------------- 86 86 INTEGER :: ji, jj, jc 87 REAL(wp) :: zQabs, zd elta, zfr87 REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus 88 88 !!--------------------------------------------------------------------- 89 89 DO jj = 1, jpj … … 91 91 92 92 zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 93 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 94 95 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 96 97 DO jc = 1, 4 ! because implicit in terms of zdelta... 98 zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) / !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 93 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt... 94 95 zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K] 96 zqlat = pQlat(ji,jj) 97 zus = pustar(ji,jj) 98 99 100 zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus ) 101 102 DO jc = 1, 4 ! because implicit in terms of zdlt... 103 zfr = MAX( 0.137_wp + 11._wp*zdlt & 104 & - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) & 105 & , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) 106 ! !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 99 107 zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 100 zd elta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj))108 zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus ) 101 109 END DO 102 110 103 dT_cs(ji,jj) = zQabs*zd elta/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)111 dT_cs(ji,jj) = zQabs*zdlt/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 104 112 105 113 END DO … … 132 140 REAL(wp) :: zdTwl, zHwl, zQabs, zfr 133 141 REAL(wp) :: zqac, ztac 134 REAL(wp) :: zal pha, zcd1, zcd2, flg142 REAL(wp) :: zalfa, zcd1, zcd2, flg 135 143 !!--------------------------------------------------------------------- 136 144 … … 142 150 !! INITIALIZATION: 143 151 zQabs = 0._wp ! total heat flux absorped in warm layer 144 zfr = zfr0 ! initial value of solar flux absorption ! LOLO: save it and use previous value !!!152 zfr = zfr0 ! initial value of solar flux absorption !#LB: save it and use previous value !!! 145 153 146 154 IF( .NOT. ln_dm2dc ) CALL sbc_dcy_param() ! we need to call sbc_dcy_param (sbcdcy.F90) because rdawn_dcy and rdusk_dcy are unkonwn otherwize! … … 148 156 ztime = REAL(nsec_day,wp)/(24._wp*3600._wp) ! time of current time step since 00:00 for current day (UTC) -> ztime = 0 -> 00:00 / ztime = 0.5 -> 12:00 ... 149 157 150 IF (lwp) THEN151 WRITE(numout,*) ''152 WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime153 END IF154 155 158 DO jj = 1, jpj 156 159 DO ji = 1, jpi … … 166 169 167 170 !***** variables for warm layer *** 168 zal pha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water(SST accurate enough!)169 170 zcd1 = SQRT(2._wp*rich*rCp0_w/(zal pha*grav*rho0_w)) !mess-o-constants 1171 zcd2 = SQRT(2._wp*zal pha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2172 173 171 zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!) 172 173 zcd1 = SQRT(2._wp*rich*rCp0_w/(zalfa*grav*rho0_w)) !mess-o-constants 1 174 zcd2 = SQRT(2._wp*zalfa*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 175 176 174 177 znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp ) ! 0<rnoon<1. => rnoon*24 = UTC time of local noon 175 178 zmidn = MOD( znoon-0.5_wp , 1._wp ) 176 zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 177 178 IF 179 zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 180 181 IF( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN 179 182 ! Dawn reset to 0! 180 183 l_exit = .TRUE. 181 184 l_destroy_wl = .TRUE. 182 END 183 184 IF 185 ENDIF 186 187 IF( .NOT. l_exit ) THEN 185 188 !! Initial test on initial guess of absorbed heat flux in warm-layer: 186 zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 187 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 188 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!! 189 190 IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 189 zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer 190 ! ! => #LB: depends of zfr, which is wild guess... Wrong!!! 191 IF( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 191 192 ! We have not started to build a WL yet (dT==0) and there's no way it can occur now 192 193 ! since zQabs <= 0._wp 193 194 ! => no need to go further 194 195 l_exit = .TRUE. 195 END 196 197 END 196 ENDIF 197 198 ENDIF 198 199 199 200 ! Okay test on updated absorbed flux: 200 ! LOLO: remove??? has a strong influence !!!201 IF ( (.NOT.(l_exit)) .AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN201 !#LB: remove??? has a strong influence !!! 202 IF( (.NOT. l_exit).AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 202 203 l_exit = .TRUE. 203 204 l_destroy_wl = .TRUE. 204 END 205 206 207 IF 205 ENDIF 206 207 208 IF( .NOT. l_exit) THEN 208 209 209 210 ! Two possibilities at this point: … … 217 218 !! some part is useless if Qsw=0 !!! 218 219 DO jl = 1, 5 219 zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 220 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 221 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) 220 zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) 222 221 zqac = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed 223 IF 222 IF( zqac <= 0._wp ) EXIT 224 223 zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 225 224 END DO 226 225 227 IF 226 IF( zqac <= 0._wp ) THEN 228 227 l_destroy_wl = .TRUE. 229 228 l_exit = .TRUE. … … 234 233 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (zdTwl = zdTwl) | 0 when gdept_1d(1)<zHwl (zdTwl = zdTwl*gdept_1d(1)/zHwl) 235 234 zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl ) 236 END 237 238 END IF !IF( .NOT. l_exit)239 240 IF 235 ENDIF 236 237 ENDIF !IF( .NOT. l_exit) 238 239 IF( l_destroy_wl ) THEN 241 240 zdTwl = 0._wp 242 241 zfr = 0.75_wp … … 244 243 zqac = 0._wp 245 244 ztac = 0._wp 246 END 247 248 IF 245 ENDIF 246 247 IF( iwait == 0 ) THEN 249 248 !! Iteration loop within bulk algo is over, time to update what needs to be updated: 250 249 dT_wl(ji,jj) = zdTwl … … 252 251 Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 253 252 Tau_ac(ji,jj) = ztac 254 END 253 ENDIF 255 254 256 255 END DO … … 276 275 REAL(wp) :: delta_skin_layer 277 276 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 278 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 277 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] 278 ! ! => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 279 279 REAL(wp), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 280 280 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] … … 282 282 REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp 283 283 !!--------------------------------------------------------------------- 284 285 zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay!284 285 zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! #LB: Double check sign + division by palpha !!! units are okay! 286 286 287 287 ztf = 0.5_wp + SIGN(0.5_wp, zQd) ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) … … 299 299 END FUNCTION delta_skin_layer 300 300 301 302 FUNCTION frac_solar_abs( pHwl ) 303 !!--------------------------------------------------------------------- 304 !! Fraction of solar heat flux absorbed inside warm layer 305 !!--------------------------------------------------------------------- 306 REAL(wp) :: frac_solar_abs 307 REAL(wp), INTENT(in) :: pHwl ! thickness of warm-layer [m] 308 !!--------------------------------------------------------------------- 309 frac_solar_abs = 1._wp - ( 0.28*0.014 *(1._wp - EXP(-pHwl/0.014)) & 310 & + 0.27*0.357*(1._wp - EXP(-pHwl/0.357)) & 311 & + 0.45*12.82*(1-EXP(-pHwl/12.82)) ) / pHwl 312 END FUNCTION frac_solar_abs 313 301 314 !!====================================================================== 302 315 END MODULE sbcblk_skin_coare -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
r11826 r11963 2 2 !!====================================================================== 3 3 !! *** MODULE sbcblk_skin_ecmwf *** 4 !! Computes: 5 !! * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer 6 !! scheme used at ECMWF 7 !! Using formulation/param. of IFS of ECMWF (cycle 40r1) 8 !! based on IFS doc (avaible online on the ECMWF's website) 9 !! 10 !! From routine turb_ecmwf maintained and developed in AeroBulk 11 !! (https://github.com/brodeau/aerobulk) 12 !! 13 !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 4 !! 5 !! Module that gathers the cool-skin and warm-layer parameterization used 6 !! by the IFS at ECMWF (recoded from scratch => 7 !! https://github.com/brodeau/aerobulk) 8 !! 9 !! Mainly based on Zeng & Beljaars, 2005 with the more recent add-up from 10 !! Takaya et al., 2010 when it comes to the warm-layer parameterization 11 !! (contribution of extra mixing due to Langmuir circulation) 12 !! 13 !! - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface skin 14 !! temperature for modeling and data assimilation. Geophysical Research 15 !! Letters, 32 (14) , pp. 1-4. 16 !! 17 !! - Takaya, Y., J.-R. Bildot, A. C. M. Beljaars, and P. A. E. M. Janssen, 18 !! 2010: Refinements to a prognostic scheme of skin sea surface 19 !! temperature. J. Geophys. Res., 115, C06009, doi:10.1029/2009JC005985 20 !! 21 !! Most of the formula are taken from the documentation of IFS of ECMWF 22 !! (cycle 40r1) (avaible online on the ECMWF's website) 23 !! 24 !! Routine 'sbcblk_skin_ecmwf' also maintained and developed in AeroBulk (as 25 !! 'mod_skin_ecmwf') 26 !! (https://github.com/brodeau/aerobulk) 27 !! 28 !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 14 29 !!---------------------------------------------------------------------- 15 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 16 !!---------------------------------------------------------------------- 17 18 !!---------------------------------------------------------------------- 19 !! cswl_ecmwf : computes the surface skin temperature (aka SSST) 20 !! based on the cool-skin/warm-layer scheme used at ECMWF 30 !! History : 4.x ! 2019-11 (L.Brodeau) Original code 21 31 !!---------------------------------------------------------------------- 22 32 USE oce ! ocean dynamics and tracers … … 25 35 USE sbc_oce ! Surface boundary condition: ocean fields 26 36 27 USE sbcblk_phy ! LOLO: all thermodynamics functions, rho_air, q_sat, etc...37 USE sbcblk_phy ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc) 28 38 29 39 USE lib_mpp ! distribued memory computing library … … 31 41 USE lib_fortran ! to use key_nosignedzero 32 42 33 34 43 IMPLICIT NONE 35 44 PRIVATE … … 38 47 39 48 !! Cool-skin related parameters: 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 41 & dT_cs !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect 50 ! ! => temperature difference between air-sea interface (z=0) 51 ! ! and right below viscous layer (z=delta) 42 52 43 53 !! Warm-layer related parameters: 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 45 & dT_wl, & !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 46 & Hz_wl !: depth of warm-layer [m] 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect 55 ! ! => difference between "almost surface (right below 56 ! ! viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 57 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m] 47 58 ! 48 59 REAL(wp), PARAMETER, PUBLIC :: rd0 = 3. !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) … … 62 73 !! Cool-skin parameterization, based on Fairall et al., 1996: 63 74 !! 64 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 65 !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 66 !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 67 !! doi:10.1029/95JC03190. 75 !! - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface 76 !! skin temperature for modeling and data assimilation. Geophysical 77 !! Research Letters, 32 (14) , pp. 1-4. 68 78 !! 69 79 !!------------------------------------------------------------------ … … 81 91 !!--------------------------------------------------------------------- 82 92 INTEGER :: ji, jj, jc 83 REAL(wp) :: zQabs, zd elta, zfr93 REAL(wp) :: zQabs, zdlt, zfr, zalfa, zus 84 94 !!--------------------------------------------------------------------- 85 95 DO jj = 1, jpj … … 87 97 88 98 zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 89 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 90 91 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 92 93 DO jc = 1, 4 ! because implicit in terms of zdelta... 94 zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 95 ! => (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 99 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt... 100 101 zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K] 102 zus = pustar(ji,jj) 103 104 zdlt = delta_skin_layer( zalfa, zQabs, zus ) 105 106 DO jc = 1, 4 ! because implicit in terms of zdlt... 107 zfr = MAX( 0.065_wp + 11._wp*zdlt & 108 & - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) & 109 & , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 110 ! ! => (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 96 111 zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 97 zd elta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj))112 zdlt = delta_skin_layer( zalfa, zQabs, zus ) 98 113 END DO 99 114 100 dT_cs(ji,jj) = zQabs*zd elta/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!)115 dT_cs(ji,jj) = zQabs*zdlt/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 101 116 102 117 END DO … … 106 121 107 122 108 109 123 SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST, pustk ) 110 124 !!--------------------------------------------------------------------- 111 125 !! 112 !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 113 !! " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 126 !! Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) with the 127 !! more recent add-up from Takaya et al., 2010 when it comes to the 128 !! warm-layer parameterization (contribution of extra mixing due to 129 !! Langmuir circulation) 130 !! 131 !! - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface skin 132 !! temperature for modeling and data assimilation. Geophysical Research 133 !! Letters, 32 (14) , pp. 1-4. 134 !! 135 !! - Takaya, Y., J.-R. Bildot, A. C. M. Beljaars, and P. A. E. M. Janssen, 136 !! 2010: Refinements to a prognostic scheme of skin sea surface 137 !! temperature. J. Geophys. Res., 115, C06009, doi:10.1029/2009JC005985 114 138 !! 115 139 !! STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER! 116 140 !! 117 !! As included in IFS Cy45r1 / E.C.M.W.F.118 141 !! ------------------------------------------------------------------ 119 142 !! … … 133 156 INTEGER :: ji, jj, jc 134 157 ! 135 REAL(wp) :: & 136 & zHwl, & !: thickness of the warm-layer [m] 137 & ztcorr, & !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 138 & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 139 & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 140 & zfr, zeta, & 141 & zusw, zusw2, & 142 & zLa, zfLa, & 143 & flg, zwf, zQabs, & 144 & ZA, ZB, zL1, zL2, & 145 & zcst0, zcst1, zcst2, zcst3 158 REAL(wp) :: zHwl !: thickness of the warm-layer [m] 159 REAL(wp) :: ztcorr !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 160 REAL(wp) :: zalfa !: thermal expansion coefficient of sea-water [1/K] 161 REAL(wp) :: zdTwl_b, zdTwl_n !: temp. diff. between "almost surface (right below viscous layer) and bottom of WL 162 REAL(wp) :: zfr, zeta 163 REAL(wp) :: zusw, zusw2 164 REAL(wp) :: zLa, zfLa 165 REAL(wp) :: flg, zwf, zQabs 166 REAL(wp) :: ZA, ZB, zL1, zL2 167 REAL(wp) :: zcst0, zcst1, zcst2, zcst3 146 168 ! 147 169 LOGICAL :: l_pustk_known … … 149 171 150 172 l_pustk_known = .FALSE. 151 IF 173 IF( PRESENT(pustk) ) l_pustk_known = .TRUE. 152 174 153 175 DO jj = 1, jpj … … 158 180 159 181 !! Previous value of dT / warm-layer, adapted to depth: 160 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) 182 flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$ 161 183 ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl 162 184 zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp ) … … 166 188 !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl ! 167 189 168 zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 169 190 zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!) 170 191 171 192 ! *** zfr = Fraction of solar radiation absorbed in warm layer (-) … … 178 199 179 200 ! Langmuir: 180 IF 201 IF( l_pustk_known ) THEN 181 202 zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6)) 182 203 ELSE 183 204 zla = 0.3_wp 184 END 205 ENDIF 185 206 zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp ) ! Eq.(6) 186 207 187 208 zwf = 0.5_wp + SIGN(0.5_wp, zQabs) ! zQabs > 0. => 1. / zQabs < 0. => 0. 188 209 189 zcst1 = vkarmn*grav*zal pha210 zcst1 = vkarmn*grav*zalfa 190 211 191 212 ! 1/L when zQabs > 0 : 192 213 zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw) 193 214 194 215 zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 ) !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ??? 195 216 196 217 zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl 197 218 198 219 ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w ) 199 220 … … 203 224 !! the sake of optimizations... So all these operations are not done 204 225 !! over and over within the iteration loop... 205 226 206 227 !! T R U L L Y I M P L I C I T => needs itteration 207 228 !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0.. … … 209 230 zdTwl_n = zdTwl_b 210 231 DO jc = 1, 10 211 232 212 233 zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence 213 234 214 235 ! 1/L when zdTwl > 0 .AND. zQabs < 0 : 215 zL1 = SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw 216 !zL1 = vkarmn*SQRT( zdTwl_n *grav*zalpha / ( 5._wp*zHwl ) ) / zusw ! => vkarmn outside, not inside zcst1 (just for this particular line) ??? 217 236 zL1 = SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalfa/( 5._wp*zHwl ) ) / zusw 237 218 238 ! Stability parameter (z/L): 219 239 zeta = (1._wp - zwf) * zHwl*zL1 + zwf * zHwl*zL2 … … 224 244 225 245 END DO 226 246 227 247 !! Update: 228 248 dT_wl(ji,jj) = zdTwl_n * ztcorr 229 249 230 250 END DO 231 251 END DO 232 252 233 253 END SUBROUTINE WL_ECMWF 234 235 254 236 255 … … 249 268 REAL(wp) :: delta_skin_layer 250 269 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 251 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 270 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] 271 ! ! => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 252 272 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] 253 273 !!--------------------------------------------------------------------- … … 255 275 !!--------------------------------------------------------------------- 256 276 ztf = 0.5_wp + SIGN(0.5_wp, pQd) ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) 257 ! ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux)258 ! 277 ! ! Qabs > 0 => warming of the viscous layer => ztf = 1 278 ! ! (ex: weak evaporation and strong positive sensible heat flux) 259 279 zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water 260 280 zusw2 = zusw*zusw … … 281 301 REAL(wp) :: ztf, zzt2 282 302 !!--------------------------------------------------------------------- 283 !284 303 zzt2 = pzeta*pzeta 285 !286 304 ztf = 0.5_wp + SIGN(0.5_wp, pzeta) ! zeta > 0 => ztf = 1 287 305 ! ! zeta < 0 => ztf = 0 288 306 PHI = ztf * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) & ! zeta > 0 289 307 & + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) ) ! zeta < 0 290 !291 308 END FUNCTION PHI 292 309 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbccpl.F90
r11831 r11963 453 453 CASE( 'conservative' ) 454 454 srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE. 455 IF 455 IF( k_ice <= 1 ) srcv(jpr_ievp)%laction = .FALSE. 456 456 CASE( 'oce and ice' ) ; srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE. 457 457 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' ) … … 558 558 srcv(jpr_botm )%clname = 'OBotMlt' 559 559 IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN 560 IF 560 IF( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN 561 561 srcv(jpr_topm:jpr_botm)%nct = nn_cats_cpl 562 562 ELSE … … 569 569 ! ! ------------------------- ! 570 570 srcv(jpr_ts_ice)%clname = 'OTsfIce' ! needed by Met Office 571 IF 572 IF 573 IF 571 IF( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' ) srcv(jpr_ts_ice)%laction = .TRUE. 572 IF( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' ) srcv(jpr_ts_ice)%nct = nn_cats_cpl 573 IF( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = nn_cats_cpl 574 574 575 575 ! ! ------------------------- ! … … 693 693 ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere 694 694 DO jn = 1, jprcv 695 IF 695 IF( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname)) 696 696 END DO 697 697 ! … … 720 720 ! =================================================== ! 721 721 DO jn = 1, jprcv 722 IF 722 IF( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) ) 723 723 END DO 724 724 ! Allocate taum part of frcv which is used even when not received as coupling field 725 IF 725 IF( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) ) 726 726 ! Allocate w10m part of frcv which is used even when not received as coupling field 727 IF 727 IF( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) ) 728 728 ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field 729 IF 730 IF 729 IF( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) ) 730 IF( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) ) 731 731 ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE. 732 732 IF( k_ice /= 0 ) THEN 733 IF 734 IF 735 END 733 IF( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) ) 734 IF( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) ) 735 ENDIF 736 736 737 737 ! ================================ ! … … 757 757 CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice' ) 758 758 ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE. 759 IF 759 IF( TRIM( sn_snd_temp%clcat ) == 'yes' ) ssnd(jps_tice)%nct = nn_cats_cpl 760 760 CASE( 'mixed oce-ice' ) ; ssnd( jps_tmix )%laction = .TRUE. 761 761 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' ) … … 777 777 ! 1. sending mixed oce-ice albedo or 778 778 ! 2. receiving mixed oce-ice solar radiation 779 IF 779 IF( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN 780 780 CALL oce_alb( zaos, zacs ) 781 781 ! Due to lack of information on nebulosity : mean clear/overcast sky … … 796 796 ssnd(jps_fice1)%laction = .TRUE. ! First-order regridded ice concentration, to be used producing atmos-to-ice fluxes (Met Office requirement) 797 797 ! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now 798 IF 799 IF 798 IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = nn_cats_cpl 799 IF( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = nn_cats_cpl 800 800 ENDIF 801 801 802 IF 802 IF(TRIM( sn_snd_ifrac%cldes ) == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 803 803 804 804 SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) … … 806 806 CASE( 'ice and snow' ) 807 807 ssnd(jps_hice:jps_hsnw)%laction = .TRUE. 808 IF 808 IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN 809 809 ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl 810 810 ENDIF 811 811 CASE ( 'weighted ice and snow' ) 812 812 ssnd(jps_hice:jps_hsnw)%laction = .TRUE. 813 IF 813 IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl 814 814 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' ) 815 815 END SELECT … … 828 828 ssnd(jps_a_p)%laction = .TRUE. 829 829 ssnd(jps_ht_p)%laction = .TRUE. 830 IF 830 IF( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 831 831 ssnd(jps_a_p)%nct = nn_cats_cpl 832 832 ssnd(jps_ht_p)%nct = nn_cats_cpl 833 833 ELSE 834 IF 834 IF( nn_cats_cpl > 1 ) THEN 835 835 CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' ) 836 836 ENDIF … … 839 839 ssnd(jps_a_p)%laction = .TRUE. 840 840 ssnd(jps_ht_p)%laction = .TRUE. 841 IF 841 IF( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 842 842 ssnd(jps_a_p)%nct = nn_cats_cpl 843 843 ssnd(jps_ht_p)%nct = nn_cats_cpl … … 914 914 CASE ( 'ice only' ) 915 915 ssnd(jps_ttilyr)%laction = .TRUE. 916 IF 916 IF( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN 917 917 ssnd(jps_ttilyr)%nct = nn_cats_cpl 918 918 ELSE 919 IF 919 IF( nn_cats_cpl > 1 ) THEN 920 920 CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_ttilyr%cldes if not exchanging category fields' ) 921 921 ENDIF … … 923 923 CASE ( 'weighted ice' ) 924 924 ssnd(jps_ttilyr)%laction = .TRUE. 925 IF 925 IF( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl 926 926 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_ttilyr%cldes;'//sn_snd_ttilyr%cldes ) 927 927 END SELECT … … 933 933 CASE ( 'ice only' ) 934 934 ssnd(jps_kice)%laction = .TRUE. 935 IF 935 IF( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN 936 936 ssnd(jps_kice)%nct = nn_cats_cpl 937 937 ELSE 938 IF 938 IF( nn_cats_cpl > 1 ) THEN 939 939 CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' ) 940 940 ENDIF … … 942 942 CASE ( 'weighted ice' ) 943 943 ssnd(jps_kice)%laction = .TRUE. 944 IF 944 IF( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl 945 945 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes;'//sn_snd_cond%cldes ) 946 946 END SELECT … … 1003 1003 ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere 1004 1004 DO jn = 1, jpsnd 1005 IF 1005 IF( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname)) 1006 1006 END DO 1007 1007 ! … … 1030 1030 CALL cpl_define(jprcv, jpsnd, nn_cplmodel) 1031 1031 1032 IF 1032 IF(ln_usecplmask) THEN 1033 1033 xcplmask(:,:,:) = 0. 1034 1034 CALL iom_open( 'cplmask', inum ) … … 1266 1266 1267 1267 IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:) ! correct this later (read from restart if possible) 1268 END 1268 ENDIF 1269 1269 ! 1270 1270 IF( ln_sdw ) THEN ! Stokes Drift correction activated … … 1415 1415 ELSE IF( srcv(jpr_qnsmix)%laction ) THEN ; zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 1416 1416 ELSE ; zqns(:,:) = 0._wp 1417 END 1417 ENDIF 1418 1418 ! update qns over the free ocean with: 1419 1419 IF( nn_components /= jp_iam_opa ) THEN … … 1687 1687 ! --- evaporation over ice (kg/m2/s) --- ! 1688 1688 DO jl=1,jpl 1689 IF 1689 IF(sn_rcv_emp%clcat == 'yes') THEN ; zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 1690 1690 ELSE ; zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 ) ; ENDIF 1691 1691 ENDDO … … 1786 1786 CASE( 'conservative' ) ! the required fields are directly provided 1787 1787 zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 1788 IF 1788 IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 1789 1789 zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl) 1790 1790 ELSE … … 1795 1795 CASE( 'oce and ice' ) ! the total flux is computed from ocean and ice fluxes 1796 1796 zqns_tot(:,:) = ziceld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1) 1797 IF 1797 IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 1798 1798 DO jl=1,jpl 1799 1799 zqns_tot(:,: ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl) … … 1897 1897 #endif 1898 1898 ! outputs 1899 IF 1900 IF 1901 IF 1902 IF 1899 IF( srcv(jpr_cal)%laction ) CALL iom_put('hflx_cal_cea' , - frcv(jpr_cal)%z3(:,:,1) * rLfus ) ! latent heat from calving 1900 IF( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * rLfus ) ! latent heat from icebergs melting 1901 IF( iom_use('hflx_rain_cea') ) CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) ) ! heat flux from rain (cell average) 1902 IF( iom_use('hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) & 1903 1903 & * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) ) ! heat flux from evap (cell average) 1904 IF 1905 IF 1904 IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) ) ! heat flux from snow (cell average) 1905 IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 1906 1906 & * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 1907 IF 1907 IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 1908 1908 & * zsnw(:,:) ) ! heat flux from snow (over ice) 1909 1909 ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp. … … 1916 1916 CASE( 'conservative' ) 1917 1917 zqsr_tot(:,: ) = frcv(jpr_qsrmix)%z3(:,:,1) 1918 IF 1918 IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 1919 1919 zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl) 1920 1920 ELSE … … 1928 1928 CASE( 'oce and ice' ) 1929 1929 zqsr_tot(:,: ) = ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1) 1930 IF 1930 IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 1931 1931 DO jl = 1, jpl 1932 1932 zqsr_tot(:,: ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl) … … 1984 1984 ! ! ========================= ! 1985 1985 CASE ('coupled') 1986 IF 1986 IF( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN 1987 1987 zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl) 1988 1988 ELSE … … 2062 2062 IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN 2063 2063 2064 IF 2064 IF( nn_components == jp_iam_opa ) THEN 2065 2065 ztmp1(:,:) = tsn(:,:,1,jp_tem) ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part 2066 2066 ELSE … … 2467 2467 IF( ssnd(jps_ficet)%laction ) THEN 2468 2468 CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 2469 END 2469 ENDIF 2470 2470 ! ! ------------------------- ! 2471 2471 ! ! Water levels to waves ! … … 2482 2482 ENDIF 2483 2483 CALL cpl_snd( jps_wlev , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 2484 END 2484 ENDIF 2485 2485 ! 2486 2486 ! Fields sent by OPA to SAS when doing OPA<->SAS coupling -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcdcy.F90
r11892 r11963 94 94 WRITE(numout,*) 95 95 ENDIF 96 END 96 ENDIF 97 97 98 98 ! Setting parameters for each new day: … … 121 121 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 122 122 ztmpm = zupusd - zlousd 123 IF 123 IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 124 124 ! 125 125 ELSE ! day time in two parts … … 135 135 ztmpm = ztmpm1 + ztmpm2 136 136 zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 137 IF 137 IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1 138 138 ENDIF 139 139 ELSE ! 24h light or 24h night … … 206 206 DO jj = 1, jpj 207 207 DO ji = 1, jpi 208 IF 208 IF( ABS(rab(ji,jj)) < 1._wp ) THEN ! day duration is less than 24h 209 209 ! When is it night? 210 210 ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 211 211 ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx ) 212 212 ! is it dawn or dusk? 213 IF 213 IF( ztest > 0._wp ) THEN 214 214 rdawn_dcy(ji,jj) = ztx 215 215 rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) ) … … 232 232 DO jj = 1, jpj 233 233 DO ji = 1, jpi 234 IF 234 IF( ABS(rab(ji,jj)) < 1._wp ) THEN ! day duration is less than 24h 235 235 rscal(ji,jj) = 0.0_wp 236 IF 236 IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN ! day time in one part 237 237 IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN 238 238 rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) … … 247 247 ENDIF 248 248 ELSE 249 IF 249 IF( raa(ji,jj) > rbb(ji,jj) ) THEN ! 24h day 250 250 rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 251 251 rscal(ji,jj) = 1._wp / rscal(ji,jj) -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcice_cice.F90
r11831 r11963 132 132 IF ( ksbc == jp_flx ) THEN 133 133 CALL cice_sbc_force(kt) 134 ELSE IF 134 ELSE IF( ksbc == jp_purecpl ) THEN 135 135 CALL sbc_cpl_ice_flx( fr_i ) 136 136 ENDIF … … 140 140 CALL cice_sbc_out ( kt, ksbc ) 141 141 142 IF 142 IF( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1) 143 143 144 144 ENDIF ! End sea-ice time step only … … 168 168 ! there is no restart file. 169 169 ! Values from a CICE restart file would overwrite this 170 IF 170 IF( .NOT. ln_rstart ) THEN 171 171 CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 172 172 ENDIF … … 177 177 178 178 ! Do some CICE consistency checks 179 IF 180 IF 179 IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 180 IF( calc_strair .OR. calc_Tsfc ) THEN 181 181 CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' ) 182 182 ENDIF 183 ELSEIF 184 IF 183 ELSEIF(ksbc == jp_blk) THEN 184 IF( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN 185 185 CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' ) 186 186 ENDIF … … 202 202 203 203 CALL cice2nemo(aice,fr_i, 'T', 1. ) 204 IF 204 IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 205 205 DO jl=1,ncat 206 206 CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. ) … … 297 297 ! forced and coupled case 298 298 299 IF 299 IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN 300 300 301 301 ztmpn(:,:,:)=0.0 … … 322 322 323 323 ! Surface downward latent heat flux (CI_5) 324 IF 324 IF(ksbc == jp_flx) THEN 325 325 DO jl=1,ncat 326 326 ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) … … 332 332 DO jj=1,jpj 333 333 DO ji=1,jpi 334 IF 334 IF(fr_i(ji,jj).eq.0.0) THEN 335 335 DO jl=1,ncat 336 336 ztmpn(ji,jj,jl)=0.0 … … 351 351 ! GBM conductive flux through ice (CI_6) 352 352 ! Convert to GBM 353 IF 353 IF(ksbc == jp_flx) THEN 354 354 ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl) 355 355 ELSE … … 360 360 ! GBM surface heat flux (CI_7) 361 361 ! Convert to GBM 362 IF 362 IF(ksbc == jp_flx) THEN 363 363 ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 364 364 ELSE … … 368 368 ENDDO 369 369 370 ELSE IF 370 ELSE IF(ksbc == jp_blk) THEN 371 371 372 372 ! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself) … … 546 546 ! Freshwater fluxes 547 547 548 IF 548 IF(ksbc == jp_flx) THEN 549 549 ! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip) 550 550 ! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below … … 552 552 ! Better to use evap and tprecip? (but for now don't read in evap in this case) 553 553 emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:)) 554 ELSE IF 554 ELSE IF(ksbc == jp_blk) THEN 555 555 emp(:,:) = (1.0-fr_i(:,:))*emp(:,:) 556 ELSE IF 556 ELSE IF(ksbc == jp_purecpl) THEN 557 557 ! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 558 558 ! This is currently as required with the coupling fields from the UM atmosphere … … 584 584 ! Scale qsr and qns according to ice fraction (bulk formulae only) 585 585 586 IF 586 IF(ksbc == jp_blk) THEN 587 587 qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:)) 588 588 qns(:,:)=qns(:,:)*(1.0-fr_i(:,:)) 589 589 ENDIF 590 590 ! Take into account snow melting except for fully coupled when already in qns_tot 591 IF 591 IF(ksbc == jp_purecpl) THEN 592 592 qsr(:,:)= qsr_tot(:,:) 593 593 qns(:,:)= qns_tot(:,:) … … 624 624 625 625 CALL cice2nemo(aice,fr_i,'T', 1. ) 626 IF 626 IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN 627 627 DO jl=1,ncat 628 628 CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. ) … … 879 879 ! B. Gather pn into global array (png) 880 880 881 IF 881 IF( jpnij > 1) THEN 882 882 CALL mppsync 883 883 CALL mppgather (pn,0,png) … … 892 892 ! (may be OK but not 100% sure) 893 893 894 IF 894 IF(nproc==0) THEN 895 895 ! pcg(:,:)=0.0 896 896 DO jn=1,jpnij … … 1015 1015 ! the lbclnk call on pn will replace these with sensible values 1016 1016 1017 IF 1017 IF(nproc==0) THEN 1018 1018 png(:,:,:)=0.0 1019 1019 DO jn=1,jpnij … … 1028 1028 ! C. Scatter png into NEMO field (pn) for each processor 1029 1029 1030 IF 1030 IF( jpnij > 1) THEN 1031 1031 CALL mppsync 1032 1032 CALL mppscatter (png,0,pn) -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcisf.F90
r11831 r11963 303 303 ! 304 304 ! Allocate public variable 305 IF 305 IF( sbc_isf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 306 306 ! 307 307 ! initialisation … … 440 440 !! Initialize arrays to 0 (each step) 441 441 zt_sum = 0.e0_wp 442 IF 442 IF( ik > 1 ) THEN 443 443 ! 1. -----------the average temperature between 200m and 600m --------------------- 444 444 DO jk = misfkt(ji,jj),misfkb(ji,jj) … … 459 459 ELSE 460 460 qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 461 END 461 ENDIF 462 462 END DO 463 463 END DO … … 496 496 ! coeficient for linearisation of potential tfreez 497 497 ! Crude approximation for pressure (but commonly used) 498 IF 498 IF( l_useCT ) THEN ! linearisation from Jourdain et al. (2017) 499 499 zlamb1 =-0.0564_wp 500 500 zlamb2 = 0.0773_wp … … 558 558 ! compute s freeze 559 559 zsfrz=(-zbqe-SQRT(zdis))*zaqer 560 IF 560 IF( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 561 561 562 562 ! compute t freeze (eq. 22) … … 578 578 579 579 ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 580 IF 580 IF( nn_gammablk < 2 ) THEN ; lit = .FALSE. 581 581 ELSE 582 582 ! check total number of iteration 583 IF 583 IF(nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 584 584 ELSE ; nit = nit + 1 585 END 585 ENDIF 586 586 587 587 ! compute error between 2 iterations 588 588 ! if needed save gammat and compute zhtflx_b for next iteration 589 589 zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 590 IF 590 IF( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 591 591 ELSE ; zhtflx_b(:,:) = zhtflx(:,:) 592 END 593 END 592 ENDIF 593 ENDIF 594 594 END DO 595 595 ! … … 718 718 pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 719 719 pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 720 END 720 ENDIF 721 721 END DO 722 722 END DO … … 757 757 ! determine the deepest level influenced by the boundary layer 758 758 DO jk = ikt+1, mbku(ji,jj) 759 IF 759 IF( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 760 760 END DO 761 761 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 789 789 ! determine the deepest level influenced by the boundary layer 790 790 DO jk = ikt+1, mbkv(ji,jj) 791 IF 791 IF( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 792 792 END DO 793 793 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 869 869 ! determine the deepest level influenced by the boundary layer 870 870 DO jk = ikt, mbkt(ji,jj) 871 IF 871 IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 872 872 END DO 873 873 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. … … 879 879 END DO 880 880 END DO 881 END 881 ENDIF 882 882 ! 883 883 !== ice shelf melting distributed over several levels ==! -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcmod.F90
r11831 r11963 243 243 fwfisf (:,:) = 0._wp ; risf_tsc (:,:,:) = 0._wp 244 244 fwfisf_b(:,:) = 0._wp ; risf_tsc_b(:,:,:) = 0._wp 245 END 245 ENDIF 246 246 IF( nn_ice == 0 ) THEN !* No sea-ice in the domain : ice fraction is always zero 247 247 IF( nn_components /= jp_iam_opa ) fr_i(:,:) = 0._wp ! except for OPA in SAS-OPA coupled case … … 399 399 emp_b (:,:) = emp (:,:) 400 400 sfx_b (:,:) = sfx (:,:) 401 IF 401 IF( ln_rnf ) THEN 402 402 rnf_b (:,: ) = rnf (:,: ) 403 403 rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) … … 437 437 IF( ln_mixcpl ) CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice ) ! forced-coupled mixed formulation after forcing 438 438 ! 439 IF 439 IF( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( ) ! Wind stress provided by waves 440 440 ! 441 441 ! !== Misc. Options ==! … … 471 471 !!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. 472 472 !!$ CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) 473 IF 473 IF( ll_wd ) THEN ! If near WAD point limit the flux for now 474 474 zthscl = atanh(rn_wd_sbcfra) ! taper frac default is .999 475 475 zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1 ! do this calc of water -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcrnf.F90
r11831 r11963 439 439 ! ! - mixed upstream-centered (ln_traadv_cen2=T) 440 440 ! 441 IF 441 IF( ln_rnf_depth ) CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already', & 442 442 & 'be spread through depth by ln_rnf_depth' ) 443 443 ! -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbctide.F90
r10068 r11963 72 72 ! Temporarily set nsec_day to beginning of day. 73 73 nsec_day_orig = nsec_day 74 IF 74 IF( nsec_day /= NINT(0.5_wp * rdt) ) THEN 75 75 kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 76 76 nsec_day = NINT(0.5_wp * rdt) -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/tideini.F90
r11831 r11963 68 68 ! 69 69 IF( ln_tide ) THEN 70 IF 70 IF(lwp) THEN 71 71 WRITE(numout,*) 72 72 WRITE(numout,*) 'tide_init : Initialization of the tidal components' … … 127 127 kt_tide = nit000 128 128 ! 129 IF 129 IF(.NOT.ln_scal_load ) rn_scal_load = 0._wp 130 130 ! 131 131 END SUBROUTINE tide_init
Note: See TracChangeset
for help on using the changeset viewer.