Changeset 5060
- Timestamp:
- Jul 16, 2024, 4:57:48 PM (8 weeks ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/dyn3d/dynredem_mod.F90
r2299 r5060 4 4 PRIVATE 5 5 PUBLIC :: dynredem_write_u, dynredem_write_v, dynredem_read_u, err 6 PUBLIC :: cre_var, get_var1,put_var1, put_var2, fil, modname, msg6 PUBLIC :: cre_var, put_var1, put_var2, fil, modname, msg 7 7 include "dimensions.h" 8 8 include "paramet.h" -
LMDZ6/trunk/libf/dyn3dmem/dynredem_mod.F90
r2299 r5060 7 7 PRIVATE 8 8 PUBLIC :: dynredem_write_u, dynredem_write_v, dynredem_read_u, err 9 PUBLIC :: cre_var, get_var1,put_var, fil, modname, msg9 PUBLIC :: cre_var, put_var, fil, modname, msg 10 10 CHARACTER(LEN=256), SAVE :: fil, modname 11 11 INTEGER, SAVE :: nvarid -
LMDZ6/trunk/libf/phylmd/cdrag_mod.F90
r4777 r5060 23 23 24 24 USE dimphy 25 USE coare_cp_mod, ONLY: coare_cp 26 USE coare30_flux_cnrm_mod, ONLY: coare30_flux_cnrm 25 27 USE indice_sol_mod 26 28 USE print_control_mod, ONLY: lunout, prt_level … … 341 343 LPWG = .false. 342 344 call ini_csts 343 call coare30_flux_cnrm(z_0m,t1(i),tsurf(i), q1(i), & 344 sqrt(zdu2),zgeop1(i)/RG,zgeop1(i)/RG,psol(i),qsurf(i),PQSAT, & 345 PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, & 346 PRESA,prain,pat1(i),z_0h, LPRECIP, LPWG, coeffs) 345 block 346 real, dimension(1) :: z0m_1d, z_0h_1d, sqrt_zdu2_1d, zgeop1_rg_1d ! convert scalar to 1D for call 347 z0m_1d = z0m 348 z_0h_1d = z0h 349 sqrt_zdu2_1d = sqrt(zdu2) 350 zgeop1_rg_1d=zgeop1(i)/RG 351 call coare30_flux_cnrm(z0m_1d,t1(i),tsurf(i), q1(i), & 352 sqrt_zdu2_1d,zgeop1_rg_1d,zgeop1_rg_1d,psol(i),qsurf(i),PQSAT, & 353 PSFTH,PFSTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI, & 354 PRESA,prain,pat1(i),z_0h_1d, LPRECIP, LPWG, coeffs) 355 356 end block 347 357 cdmm(i) = coeffs(1) 348 358 cdhh(i) = coeffs(2) -
LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90
r4722 r5060 5 5 ! ######### 6 6 7 !------------------------ Modif Olivier Torres pour rajout fonction et pas appel ---------------- 8 9 10 11 real function PSIFCTT(zet) 12 13 real, intent(in) :: zet 14 15 if(zet<0) then 16 x=(1.-(15*zet))**.5 17 psik=2*log((1+x)/2) 18 x=(1.-(34.15*zet))**.3333 19 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 20 f=zet*zet/(1+zet*zet) 21 PSIFCTT=(1-f)*psik+f*psic 22 23 else 24 c=min(50.,.35*zet) 25 PSIFCTT=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) 26 endif 27 end FUNCTION PSIFCTT 28 29 real function PSIFCTU(zet) 30 31 real, intent(in) :: zet 32 33 if (zet<0) then 34 35 x=(1.-15.*zet)**.25 36 psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) 37 x=(1.-10.15*zet)**.3333 38 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 39 f=zet*zet/(1+zet*zet) 40 PSIFCTU=(1-f)*psik+f*psic 41 else 42 c=min(50.,.35*zet) 43 PSIFCTU=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) 44 endif 45 END FUNCTION PSIFCTU 46 47 !----------------------------------- Fin Modif --------------------------------------------------- 7 module coare30_flux_cnrm_mod 8 implicit none 9 private 10 public COARE30_FLUX_CNRM 11 12 contains 48 13 49 14 50 15 SUBROUTINE COARE30_FLUX_CNRM(PZ0SEA,PTA,PSST,PQA, & 51 16 PVMOD,PZREF,PUREF,PPS,PQSATA,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,& 52 PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs ) 17 PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs ) 53 18 ! ####################################################################### 54 19 ! 55 20 ! 56 !!**** *COARE25_FLUX* 21 !!**** *COARE25_FLUX* 57 22 !! 58 23 !! PURPOSE 59 24 !! ------- 60 25 ! Calculate the surface fluxes of heat, moisture, and momentum over 61 ! sea surface with bulk algorithm COARE3.0. 62 ! 26 ! sea surface with bulk algorithm COARE3.0. 27 ! 63 28 !!** METHOD 64 29 !! ------ 65 30 ! transfer coefficients were obtained using a dataset which combined COARE 66 31 ! data with those from three other ETL field experiments, and reanalysis of 67 ! the HEXMAX data (DeCosmos et al. 1996). 68 ! ITERMAX=3 69 ! Take account of the surface gravity waves on the velocity roughness and 32 ! the HEXMAX data (DeCosmos et al. 1996). 33 ! ITERMAX=3 34 ! Take account of the surface gravity waves on the velocity roughness and 70 35 ! hence the momentum transfer coefficient 71 36 ! NGRVWAVES=0 no gravity waves action (Charnock) !default value … … 77 42 !! 78 43 !! IMPLICIT ARGUMENTS 79 !! ------------------ 80 !! 44 !! ------------------ 45 !! 81 46 !! REFERENCE 82 47 !! --------- … … 85 50 !! Gosnell et al (1995), JGR, 437-442 86 51 !! Fairall et al (1996), JGR, 1295-1308 87 !! 52 !! 88 53 !! AUTHOR 89 54 !! ------ … … 95 60 !! B. Decharme 06/2009 limitation of Ri 96 61 !! B. Decharme 09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90 97 !! B. Decharme 06/2013 bug in z0 (output) computation 62 !! B. Decharme 06/2013 bug in z0 (output) computation 98 63 !! J.Escobar 06/2013 for REAL4/8 add EPSILON management 99 64 !! C. Lebeaupin 03/2014 bug if PTA=PSST and PEXNA=PEXNS: set a minimum value … … 111 76 USE dimphy 112 77 USE indice_sol_mod 78 USE coare_cp_mod, ONLY: PSIFCTT=>psit_30, PSIFCTU=>psiuo 113 79 114 80 !-------------------------------- … … 159 125 ! 160 126 REAL, DIMENSION(klon), INTENT(INOUT) :: PZ0SEA! roughness length over the ocean 161 ! 127 ! 162 128 ! surface fluxes : latent heat, sensible heat, friction fluxes 163 129 REAL, DIMENSION(klon), INTENT(OUT) :: PSFTH ! heat flux (W/m2) … … 176 142 177 143 LOGICAL, INTENT(IN) :: LPRECIP ! 178 LOGICAL, INTENT(IN) :: LPWG ! 144 LOGICAL, INTENT(IN) :: LPWG ! 179 145 real, dimension(3), intent(inout) :: coeffs 180 146 ! … … 195 161 REAL, DIMENSION(SIZE(PTA)) :: PEXNS ! Exner function at atm level 196 162 ! 197 REAL, DIMENSION(SIZE(PTA)) :: ZO ! rougness length ref 163 REAL, DIMENSION(SIZE(PTA)) :: ZO ! rougness length ref 198 164 REAL, DIMENSION(SIZE(PTA)) :: ZWG ! gustiness factor (m/s) 199 165 ! … … 204 170 REAL, DIMENSION(SIZE(PTA)) :: ZQSR !humidity scaling parameter "qstar" (kg/kg) 205 171 ! 206 REAL, DIMENSION(SIZE(PTA)) :: ZU10,ZT10 !vertical profils (10-m height) 172 REAL, DIMENSION(SIZE(PTA)) :: ZU10,ZT10 !vertical profils (10-m height) 207 173 REAL, DIMENSION(SIZE(PTA)) :: ZVISA !kinematic viscosity of dry air 208 174 REAL, DIMENSION(SIZE(PTA)) :: ZO10,ZOT10 !roughness length at 10m … … 215 181 REAL, DIMENSION(SIZE(PTA)) :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact 216 182 ! 217 REAL, DIMENSION(SIZE(PTA)) :: ZZL,ZZTL!,ZZQL !Obukhovs stability 183 REAL, DIMENSION(SIZE(PTA)) :: ZZL,ZZTL!,ZZQL !Obukhovs stability 218 184 !param. z/l for u,T,q 219 185 REAL, DIMENSION(SIZE(PTA)) :: ZRR 220 186 REAL, DIMENSION(SIZE(PTA)) :: ZOT,ZOQ !rougness length ref 221 REAL, DIMENSION(SIZE(PTA)) :: ZPUZ,ZPTZ,ZPQZ !PHI funct. for u,T,q 187 REAL, DIMENSION(SIZE(PTA)) :: ZPUZ,ZPTZ,ZPQZ !PHI funct. for u,T,q 222 188 ! 223 189 REAL, DIMENSION(SIZE(PTA)) :: ZBF !constants to compute gustiness factor … … 235 201 REAL, DIMENSION(SIZE(PTA)) :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact 236 202 REAL, DIMENSION(SIZE(PTA)) :: ZXLR ! vaporisation heat at a given temperature 237 REAL, DIMENSION(SIZE(PTA)) :: ZCPLW ! specific heat for water at a given temperature 203 REAL, DIMENSION(SIZE(PTA)) :: ZCPLW ! specific heat for water at a given temperature 238 204 ! 239 205 REAL, DIMENSION(SIZE(PTA)) :: ZUSTAR2 ! square of friction velocity … … 254 220 REAL :: QSAT_SEAWATER 255 221 REAL :: QSATSEAW_1D 256 REAL, EXTERNAL :: PSIFCTU, PSIFCTT257 222 ! 258 223 INTEGER :: J, JLOOP !loop indice … … 264 229 265 230 REAL :: XVCHRNK = 0.021 266 REAL :: XVZ0CM = 1.0E-5 231 REAL :: XVZ0CM = 1.0E-5 267 232 !REAL :: XRIMAX 268 233 … … 320 285 ! 321 286 !------------------------------------------------------------------------------- 322 ! 2. INITIAL GUESS FOR THE ITERATIVE METHOD 287 ! 2. INITIAL GUESS FOR THE ITERATIVE METHOD 323 288 ! ------------------------------------- 324 289 ! 325 ! 2.0 Temperature 290 ! 2.0 Temperature 326 291 ! 327 292 ! Set a non-zero value for the temperature gradient 328 293 ! 329 WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.) 294 WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.) 330 295 ZTA(:)=PTA(:)-1E-3 331 296 ELSEWHERE 332 ZTA(:)=PTA(:) 297 ZTA(:)=PTA(:) 333 298 ENDWHERE 334 299 335 ! 2.1 Wind and humidity 336 ! 337 ! Sea surface specific humidity 338 ! 339 !PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:)) 340 PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:)) 341 342 343 344 ! 345 ! Set a minimum value to wind 300 ! 2.1 Wind and humidity 301 ! 302 ! Sea surface specific humidity 303 ! 304 !PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:)) 305 PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:)) 306 307 308 309 ! 310 ! Set a minimum value to wind 346 311 ! 347 312 !ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:)) … … 352 317 353 318 ! 354 ! Specific humidity at saturation at the atm. level 319 ! Specific humidity at saturation at the atm. level 355 320 ! 356 321 ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD)) 357 !ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:)) 358 ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:)) 322 !ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:)) 323 ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:)) 359 324 360 325 … … 365 330 IF (LPWG) ZWG(:) = 0.5 366 331 ! 367 ZCHARN(:) = 0.011 332 ZCHARN(:) = 0.011 368 333 ! 369 334 DO J=1,SIZE(PTA) 370 335 ! 371 336 ! 2.2 initial guess 372 ! 337 ! 373 338 ZDU(J) = ZVMOD(J) !wind speed difference with surface current(=0) (m/s) 374 339 !initial guess for gustiness factor … … 384 349 ZVISA(J) = 1.326E-5*(1.+6.542E-3*(ZTA(J)-XTT)+& 385 350 8.301E-6*(ZTA(J)-XTT)**2-4.84E-9*(ZTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11 386 ! 351 ! 387 352 ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J) 388 353 ZCD(J) = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2 !drag coefficient … … 400 365 ! &((ZTA(J)-XTT)*ZDUWG(J)**2) 401 366 ZRIBU(J) = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*ZTA(J)*ZDQ(J))/& 402 (ZTA(J)*ZDUWG(J)**2) 367 (ZTA(J)*ZDUWG(J)**2) 403 368 ! 404 369 IF (ZRIBU(J)<0.) THEN 405 370 ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J)) !Unstable G and F 406 371 ELSE 407 ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable 372 ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable 408 373 ENDIF 409 374 ! … … 413 378 ! 414 379 ! First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL) 415 ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF /ZL10))416 ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF /ZL10))417 ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF /ZL10))380 ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF(1)/ZL10(1))) 381 ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1))) 382 ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1))) 418 383 ! 419 384 ZZL(:) = 0.0 … … 431 396 IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018 432 397 ! 433 ! 3. ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 398 ! 3. ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 434 399 ! ------------------------------------------- 435 400 ! … … 441 406 ENDDO 442 407 ! 443 408 444 409 ! 445 410 DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop … … 447 412 DO J=1,SIZE(PTA) 448 413 ! 449 IF (JLOOP .GT.ITERMAX(J)) CYCLE414 IF (JLOOP>ITERMAX(J)) CYCLE 450 415 ! 451 416 IF (NGRVWAVES==0) THEN … … 453 418 ELSE IF (NGRVWAVES==1) THEN 454 419 ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 & 455 + 0.11*ZVISA(J)/ZUSR(J) !Oost et al. 2002 420 + 0.11*ZVISA(J)/ZUSR(J) !Oost et al. 2002 456 421 ELSE IF (NGRVWAVES==2) THEN 457 422 ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 & 458 + 0.11*ZVISA(J)/ZUSR(J) !Taulor and Yelland 2001 423 + 0.11*ZVISA(J)/ZUSR(J) !Taulor and Yelland 2001 459 424 ENDIF 460 425 ! … … 465 430 ZZL(J) = XKARMAN * XG * PUREF(J) * & 466 431 ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*ZTA(J)*ZQSR(J) ) / & 467 ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) ) 468 ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J) ! for T 432 ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) ) 433 ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J) ! for T 469 434 ! ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J) ! for Q 470 435 ENDDO 471 436 ! 472 ZPUZ(:) = PSIFCTU(ZZL( :))473 ZPTZ(:) = PSIFCTT(ZZTL( :))437 ZPUZ(:) = PSIFCTU(ZZL(1)) 438 ZPTZ(:) = PSIFCTT(ZZTL(1)) 474 439 ! 475 440 DO J=1,SIZE(PTA) 476 441 ! 477 ! ZPQZ(J)=PSIFCTT(ZZQL(J)) 442 ! ZPQZ(J)=PSIFCTT(ZZQL(J)) 478 443 ZPQZ(J) = ZPTZ(J) 479 444 ! … … 493 458 ZWG(J) = 0.2 494 459 ENDIF 495 ENDIF 460 ENDIF 496 461 ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2) 497 462 ! … … 515 480 ! 516 481 ! 517 ! 4. transfert coefficients PCD, PCH and PCE 518 ! and neutral PCDN, ZCHN, ZCEN 482 ! 4. transfert coefficients PCD, PCH and PCE 483 ! and neutral PCDN, ZCHN, ZCEN 519 484 ! 520 485 PCD(J) = (ZUSR(J)/ZDUWG(J))**2. … … 528 493 ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT) 529 494 ! 530 ! 4. 2 surface fluxes 531 ! 532 IF (ABS(PCDN(J))>1.E-2) THEN !!!! secure COARE3.0 CODE 495 ! 4. 2 surface fluxes 496 ! 497 IF (ABS(PCDN(J))>1.E-2) THEN !!!! secure COARE3.0 CODE 533 498 write(*,*) 'pb PCDN in COARE30: ',PCDN(J) 534 499 write(*,*) 'point: ',J,"/",SIZE(PTA) … … 543 508 ZHF(J) = PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J) 544 509 ZEF(J) = PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J) 545 ! 510 ! 546 511 ! 4.3 Contributions to surface fluxes due to rainfall 547 512 ! … … 550 515 ! sec) au cas humide. 551 516 IF (LPRECIP) THEN 552 ! 517 ! 553 518 ! heat surface fluxes 554 519 ! … … 561 526 ! 562 527 ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (ZTA(J)/XTT)**1.94 ! water vapour diffusivity from eq (13.3) 563 ! ! of Pruppacher and Klett (1978) 528 ! ! of Pruppacher and Klett (1978) 564 529 ZALFAC(J)= 1.0 / (1.0 + & ! Eq.11 in GoF95 565 ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD)) ! ZALFAC=wet-bulb factor (sans dim) 530 ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD)) ! ZALFAC=wet-bulb factor (sans dim) 566 531 ZCPLW(J) = 4224.8482 + ZTAC(J) * & 567 532 ( -4.707 + ZTAC(J) * & 568 533 (0.08499 + ZTAC(J) * & 569 534 (1.2826e-3 + ZTAC(J) * & 570 (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat 571 ! 535 (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat 536 ! 572 537 ZRF(J) = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * & !Eq.12 in GoF95 !SIGNE? 573 538 (PSST(J) - ZTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD ) 574 539 ! 575 ! Momentum flux due to rainfall 540 ! Momentum flux due to rainfall 576 541 ! 577 542 ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96 … … 580 545 ! 581 546 ! 4.4 Webb correction to latent heat flux 582 ! 547 ! 583 548 ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) & 584 - ZUSR(J)*ZTSR(J)/ZTA(J) ! Eq.21*rhoa in FBR96 585 ! 586 ! 4.5 friction velocity which contains correction du to rain 549 - ZUSR(J)*ZTSR(J)/ZTA(J) ! Eq.21*rhoa in FBR96 550 ! 551 ! 4.5 friction velocity which contains correction du to rain 587 552 ! 588 553 ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J) … … 590 555 ! 591 556 ! 4.6 Total surface fluxes 592 ! 557 ! 593 558 PSFTH (J) = ZHF(J) + ZRF(J) 594 559 PSFTQ (J) = ZEF(J) / ZLV(J) 595 ! 560 ! 596 561 ENDIF 597 ENDDO 562 ENDDO 598 563 599 564 … … 601 566 PCE,& 602 567 PCH] 603 568 604 569 !------------------------------------------------------------------------------- 605 570 ! 606 ! 5. FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS 571 ! 5. FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS 607 572 ! ----------- 608 573 ! 5.1 Richardson number 609 ! 574 ! 610 575 ! 611 576 !------------STOP LA -------------------- 612 577 !ZDIRCOSZW(:) = 1. 613 578 ! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,ZTA,ZQASAT,& 614 ! PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI ) 579 ! PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI ) 615 580 !! 616 581 !! 5.2 Aerodynamical conductance and resistance 617 !! 582 !! 618 583 !ZAC(:) = PCH(:)*ZVMOD(:) 619 584 !PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON) … … 630 595 ! 631 596 END SUBROUTINE COARE30_FLUX_CNRM 597 598 end module coare30_flux_cnrm_mod -
LMDZ6/trunk/libf/phylmd/coare_cp.F90
r4722 r5060 1 real function psit_30(zet) 2 3 real, intent(in) :: zet 4 5 if(zet<0) then 6 x=(1.-(15*zet))**.5 7 psik=2*log((1+x)/2) 8 x=(1.-(34.15*zet))**.3333 9 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 10 f=zet*zet/(1+zet*zet) 11 psit_30=(1-f)*psik+f*psic 12 13 else 14 c=min(50.,.35*zet) 15 psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) 1 module coare_cp_mod 2 implicit none 3 private 4 public psit_30, psiuo, coare_cp 5 6 contains 7 8 real function psit_30(zet) 9 implicit none 10 real, intent(in) :: zet 11 12 real :: x, psik, psic, f, c 13 14 if(zet<0) then 15 x=(1.-(15.*zet))**.5 16 psik=2.*log((1.+x)/2) 17 x=(1.-(34.15*zet))**.3333 18 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 19 f=zet*zet/(1.+zet*zet) 20 psit_30=(1.-f)*psik+f*psic 21 else 22 c=min(50.,.35*zet) 23 psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525) 24 endif 25 26 end function psit_30 27 28 real function psiuo(zet) 29 implicit none 30 real, intent(in) :: zet 31 32 real :: x, psik, psic, f, c 33 34 if (zet<0) then 35 x=(1.-15.*zet)**.25 36 psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) 37 x=(1.-10.15*zet)**.3333 38 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 39 f=zet*zet/(1+zet*zet) 40 psiuo=(1-f)*psik+f*psic 41 else 42 c=min(50.,.35*zet) 43 psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) 16 44 endif 17 end FUNCTION psit_30 18 19 real function psiuo(zet) 20 21 real, intent(in) :: zet 22 23 if (zet<0) then 24 25 x=(1.-15.*zet)**.25 26 psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.) 27 x=(1.-10.15*zet)**.3333 28 psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.) 29 f=zet*zet/(1+zet*zet) 30 psiuo=(1-f)*psik+f*psic 31 else 32 c=min(50.,.35*zet) 33 psiuo=-((1+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525) 34 endif 35 END FUNCTION psiuo 36 37 38 real function gen_q1(ts,p,dq) 39 40 real, intent(in) :: ts,p,dq 41 real es,qs 42 43 es = 6.112*exp(17.502*(ts-273.15)/(ts-32.18))*.98*(1.0007+3.46e-8*p) 44 qs = es*62.197/(p-37.8*es) 45 46 q1 = dq + qs 47 48 end function gen_q1 49 50 45 46 end function psiuo 51 47 52 48 … … 61 57 !version with shortened iteration modified Rt and Rq 62 58 !uses wave information wave period in s and wave ht in m 63 !no wave, standard coare 2.6 charnock: jwave=0 59 !no wave, standard coare 2.6 charnock: jwave=0 64 60 !Oost et al. zo=50/2/pi L (u*/c)**4.5 if jwave=1 65 61 !taylor and yelland zo=1200 h*(L/h)**4.5 jwave=2 … … 104 100 real old_usr, old_tsr, old_qsr,tmp 105 101 106 real, external :: psit_30, psiuo,grv102 real, external :: grv 107 103 integer i,j,k 108 104 !---------------- Rajout pour prendre en compte différent Z0 --------------------------------! … … 113 109 !--------------------------------------------------------------------------------------------- 114 110 115 Ribcu=-zu/zi/.004/Beta**3 116 cpa=1004.67 111 Ribcu=-zu/zi/.004/Beta**3 112 cpa=1004.67 117 113 118 114 t_c = t - 273.15 119 115 120 116 121 Le=(2.501-.00237*(t_c-dt))*1e6 122 123 124 visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) 125 126 127 cpv=cpa*(1+0.84*Q) 128 rhoa=P/(Rgas*t*(1+0.61*Q)) 129 130 131 132 133 134 135 ug=.2 136 137 ut=sqrt(du*du+ug*ug) 117 Le=(2.501-.00237*(t_c-dt))*1e6 118 119 120 visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3) 121 122 123 cpv=cpa*(1+0.84*Q) 124 rhoa=P/(Rgas*t*(1+0.61*Q)) 125 126 127 128 129 130 131 ug=.2 132 133 ut=sqrt(du*du+ug*ug) 138 134 ut=MAX(ut , 0.1 * MIN(10.,zu) ) 139 135 140 u10=ut*log(10/1e-4)/log(zu/1e-4) 136 u10=ut*log(10/1e-4)/log(zu/1e-4) 141 137 usr=.035*u10 ! turbulent friction velocity (m/s), including gustiness 142 138 zom10=0.011*usr*usr/grav+0.11*visa/usr ! roughness length for u (smith 88) 143 Cd10=(von/log(10/zom10))**2 139 Cd10=(von/log(10/zom10))**2 144 140 ! zoh10=0.40*visa/usr+1.4e-5 145 ! Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca 146 Ch10=0.00115 147 Ct10=Ch10/sqrt(Cd10) 141 ! Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca 142 Ch10=0.00115 143 Ct10=Ch10/sqrt(Cd10) 148 144 zot10=10/exp(von/Ct10) ! roughness length for t 149 Cd=(von/log(zu/zom10))**2 150 Ct=von/log(zt/zot10) 151 CC=von*Ct/Cd 145 Cd=(von/log(zu/zom10))**2 146 Ct=von/log(zt/zot10) 147 CC=von*Ct/Cd 152 148 153 149 ut0 = ut … … 155 151 156 152 ut = ut0 157 Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 158 159 if (Ribu .LT. 0) then160 zetu=CC*Ribu/(1+Ribu/Ribcu) 161 else 162 zetu=CC*Ribu*(1+27/9*Ribu/CC) 153 Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2 154 155 if (Ribu < 0) then 156 zetu=CC*Ribu/(1+Ribu/Ribcu) 157 else 158 zetu=CC*Ribu*(1+27/9*Ribu/CC) 163 159 endif 164 160 165 L10=zu/zetu 166 167 ! if (zetu .GT. 50) then 168 ! nits=1 161 L10=zu/zetu 162 163 ! if (zetu .GT. 50) then 164 ! nits=1 169 165 ! endif 170 166 171 167 usr=ut*von/(log(zu/zom10)-psiuo(zetu)) 172 tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) 173 qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) 168 tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10)) 169 qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10)) 174 170 ! tkt=.001 175 171 176 172 ! charnock constant - lin par morceau - constant 177 173 if ( ut <= 10. ) then 178 charn=0.011 174 charn=0.011 179 175 else 180 if (ut .GT.18) then176 if (ut > 18) then 181 177 charn=0.018 182 178 else 183 charn=0.011+(ut-10)/(18-10)*(0.018-0.011) 179 charn=0.011+(ut-10)/(18-10)*(0.018-0.011) 184 180 endif 185 181 endif … … 192 188 193 189 !*************** bulk loop ************ 194 do i=1, nits 190 do i=1, nits 195 191 196 192 zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q) … … 200 196 ! ELSE IF (NGRVWAVES==1) THEN 201 197 ! zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 & 202 ! + 0.11*visa/usr !Oost et al. 2002 198 ! + 0.11*visa/usr !Oost et al. 2002 203 199 ! ELSE IF (NGRVWAVES==2) THEN 204 200 ! zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 & 205 ! + 0.11*visa/usr !Taulor and Yelland 2001 206 ! ENDIF 207 208 zom=charn*usr*usr/grav+0.11*visa/usr 209 210 rr=zom*usr/visa 211 L=zu/zet 201 ! + 0.11*visa/usr !Taulor and Yelland 2001 202 ! ENDIF 203 204 zom=charn*usr*usr/grav+0.11*visa/usr 205 206 rr=zom*usr/visa 207 L=zu/zet 212 208 zoq=min(1.15e-4,5.5e-5/rr**.6) ! a modifier 213 209 zot=zoq ! a modifier 214 ! zot=0.40*visa/usr+1.4e-5 210 ! zot=0.40*visa/usr+1.4e-5 215 211 216 212 old_usr = usr … … 218 214 old_qsr = qsr 219 215 220 usr=ut*von/(log(zu/zom)-psiuo(zu/L)) 221 tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) 222 qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) 223 224 Bf=-grav/t*usr*(tsr+.61*t*qsr) 225 226 if (Bf .GT. 0) then227 ug=Beta*(Bf*zi)**.333 228 else 229 ug=.2 216 usr=ut*von/(log(zu/zom)-psiuo(zu/L)) 217 tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L)) 218 qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L)) 219 220 Bf=-grav/t*usr*(tsr+.61*t*qsr) 221 222 if (Bf > 0) then 223 ug=Beta*(Bf*zi)**.333 224 else 225 ug=.2 230 226 endif 231 227 232 ut=sqrt(du*du+ug*ug) 228 ut=sqrt(du*du+ug*ug) 233 229 ut=MAX(ut , 0.1 * MIN(10.,zu) ) 234 230 … … 237 233 238 234 ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress 239 ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr 240 ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr 235 ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr 236 ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr 241 237 242 238 tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du … … 249 245 rugosh = zot 250 246 251 252 247 end subroutine coare_cp 248 249 end module coare_cp_mod
Note: See TracChangeset
for help on using the changeset viewer.