Changeset 633 for trunk/NEMO/OPA_SRC
- Timestamp:
- 2007-03-02T17:46:40+01:00 (17 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/SBC/flx_core.h90 (modified) (42 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/SBC/flx_core.h90
r606 r633 22 22 CHARACTER (LEN=32), DIMENSION (jpfile) :: clvarname 23 23 CHARACTER (LEN=50), DIMENSION (jpfile) :: clname 24 CHARACTER (LEN=32) :: clvarkatax , clvarkatay ! variable name for katabatic masks25 CHARACTER (LEN=256) :: clnamekata ! filename for katabatic masks26 24 27 25 INTEGER :: isnap … … 32 30 freqh ! Frequency for each forcing file (hours) 33 31 ! ! a negative value of -12 corresponds to monthly 34 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! used for modif wind stress in the first sea points35 & rmskkatax, rmskkatay !: mask ocean to increase wind stress in first sea points36 32 REAL(wp), DIMENSION(jpi,jpj,jpfile,2) :: & 37 33 flxdta !: set of NCAR 6hourly/daily/monthly fluxes 38 34 LOGICAL :: & 39 & ln_2m = .FALSE. , &!: logical flag for height of air temp. and hum.40 & ln_kata = .FALSE. !: logical flag for Katabatic winds enhancement.35 & ln_2m = .FALSE. !: logical flag for height of air temp. and hum. 36 REAL(wp) :: alpha_precip=1. !: multiplication factor for precipitation 41 37 !!---------------------------------------------------------------------- 42 38 !! History : … … 48 44 !! ! - Implement reading of 6-hourly fields 49 45 !! ! 06-02 (S. Masson, G. Madec) IOM read + NEMO "style" 50 !! ! 07-06 (P. Mathiot, J-M Molines) Katabatic winds enhancement51 46 !! ! 12-06 (L. Brodeau) handle both 2m and 10m input fields 52 47 !!---------------------------------------------------------------------- … … 103 98 !! * modules used 104 99 USE iom 100 USE restart 105 101 USE par_oce 106 102 USE flx_oce … … 111 107 USE lbclnk 112 108 USE dtatem ! FOR Ce = F(SST(levitus)): 113 109 114 110 !! * arguments 115 111 INTEGER, INTENT( in ) :: kt ! ocean time step … … 128 124 zadatrjb, & ! date (noninteger day) at which we read the forcings 129 125 zxy , & ! scalar for temporal interpolation 130 zcof 126 zcof , zzu , zzv ! scalar 131 127 REAL(wp), DIMENSION(jpi,jpj, jpfile) :: & 132 128 flxnow ! input flux values at current time step … … 161 157 catm1 ! 162 158 163 NAMELIST /namcore/ clname, clvarname, freqh, ln_2m, ln_kata159 NAMELIST /namcore/ clname, clvarname, freqh, ln_2m, alpha_precip 164 160 !!--------------------------------------------------------------------- 165 161 … … 169 165 170 166 !--- calculation for monthly data 167 ! i15 = 0 if first half of current month 168 ! i15 = 1 if second half of current month 171 169 i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 172 170 iman = 12 … … 181 179 zadatrjb = adatrj - rdt / rday 182 180 ihour = INT( ( zadatrjb - INT(zadatrjb) ) * 24 ) 183 181 184 182 ! ! ------------------------ ! 185 183 IF( kt == nit000 ) THEN ! first call kt=nit000 ! 186 184 ! ! ------------------------ ! 185 CALL core_rst ( kt, 'READ' ) 187 186 188 187 !--- Initializes default values of file names, frequency of forcings 189 188 ! and variable to be read in the files, before reading namelist 190 clname(1) = 'precip _core.nc' ; freqh(1) = -12 ; clvarname(1) = 'precip' ! monthly191 clname(2) = 'u10 _core.nc' ; freqh(2) = 24 ; clvarname(2) = 'u_10'! 6 hourly192 clname(3) = 'v10 _core.nc' ; freqh(3) = 24 ; clvarname(3) = 'v_10'! 6 hourly193 clname(4) = 'q10 _core.nc' ; freqh(4) = 24 ; clvarname(4) = 'q_10'! 6 hourly194 clname(5) = ' tot_solar_core.nc'; freqh(5) = 24 ; clvarname(5) = 'tot_solar'! daily195 clname(6) = ' therm_rad_core.nc'; freqh(6) = 24 ; clvarname(6) = 'therm_rad'! daily196 clname(7) = 't emp_10m_core.nc' ; freqh(7) = 24 ; clvarname(7) = 't_10'! 6 hourly197 clname(8) = 'snow _core.nc' ; freqh(8) = -12 ; clvarname(8) = 'snow' ! monthly189 clname(1) = 'precip.nc' ; freqh(1) = -12 ; clvarname(1) = 'precip' ! monthly 190 clname(2) = 'u10.nc' ; freqh(2) = 6 ; clvarname(2) = 'u10' ! 6 hourly 191 clname(3) = 'v10.nc' ; freqh(3) = 6 ; clvarname(3) = 'v10' ! 6 hourly 192 clname(4) = 'q10.nc' ; freqh(4) = 6 ; clvarname(4) = 'q10' ! 6 hourly 193 clname(5) = 'radsw.nc' ; freqh(5) = 24 ; clvarname(5) = 'radsw' ! daily 194 clname(6) = 'radlw.nc' ; freqh(6) = 24 ; clvarname(6) = 'radlw' ! daily 195 clname(7) = 't10.nc' ; freqh(7) = 6 ; clvarname(7) = 't10' ! 6 hourly 196 clname(8) = 'snow.nc' ; freqh(8) = -12 ; clvarname(8) = 'snow' ! monthly 198 197 199 198 !--- Read Namelist namcore : OMIP/CORE 200 199 REWIND ( numnam ) 201 200 READ ( numnam, namcore ) 201 202 ! in case of ln_2m, air temp. and humidity are given at 2 m, thus file name and variable name are changed 203 IF ( ln_2m ) THEN 204 clname(4) = 'q2.nc' ; clvarname(4) = 'q2' 205 clname(7) = 't2.nc' ; clvarname(7) = 't2' 206 ! but if for some reason, clname and varname were also in the name list, we screw them up ! 207 ! re-read the namelist, to restore clname, varname to their namelist value 208 REWIND ( numnam ) 209 READ ( numnam, namcore ) 210 ENDIF 211 202 212 IF(lwp) THEN 203 213 WRITE(numout,*)' ' 204 WRITE(numout,*)' flx mod/flx_core : global CORE fields in NetCDF format'205 WRITE(numout,*)' ~~~~~~~~~~ ~~~~~'206 WRITE(numout,*) ' ln_2m = ', ln_2m207 WRITE(numout,*) ' ln_kata = ', ln_kata208 WRITE(numout,*) ' list of files and frequency (hour), or monthly (-12) '214 WRITE(numout,*)' flx_core : global CORE fields in NetCDF format' 215 WRITE(numout,*)' ~~~~~~~~~~ ' 216 WRITE(numout,*) ' ln_2m = ', ln_2m 217 WRITE(numout,*) ' alpha_precip = ', alpha_precip 218 WRITE(numout,*) ' list of files and frequency (hour), or monthly (-12) ' 209 219 DO ji = 1, jpfile 210 220 WRITE(numout,*) trim(clname(ji)), ' frequency:', freqh(ji) … … 235 245 nrecflx (:) = 0 ! switch for reading flux data for each file 236 246 nrecflx2(:) = 0 ! switch for reading flux data for each file 237 247 238 248 !--- Open the files of the list 239 249 DO ji = 1, jpfile … … 251 261 DO ji = 1, jpfile 252 262 253 IF ( freqh(ji) .GT.0 .AND. freqh(ji) .LT.24 ) THEN !--- Case of 6-hourly flux data263 IF ( freqh(ji) > 0 .AND. freqh(ji) < 24 ) THEN !--- Case of 6-hourly flux data 254 264 !--- calculate current snapshot from hour of day 255 265 isnap = ihour / INT( freqh(ji) ) + 1 … … 263 273 ENDIF 264 274 265 ELSE IF ( freqh(ji) .EQ.24 ) THEN !--- Case of daily flux data275 ELSE IF ( freqh(ji) == 24 ) THEN !--- Case of daily flux data 266 276 267 277 !--- reading is necessary when nrecflx(ji) differs from nday … … 274 284 ENDIF 275 285 276 ELSE IF ( freqh(ji) .EQ.-12 ) THEN !--- Case monthly data from CORE286 ELSE IF ( freqh(ji) == -12 ) THEN !--- Case monthly data from CORE 277 287 ! we read two months all the time although we 278 288 ! could only read one and swap the arrays … … 295 305 ENDIF 296 306 END DO ! end of loop over forcing files to verify if reading is necessary 297 307 298 308 !--- Printout if required 299 309 ! IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN … … 320 330 ! 321 331 DO ji = 1, jpfile 322 IF( freqh(ji) .EQ.-12 ) THEN332 IF( freqh(ji) == -12 ) THEN 323 333 zxy = REAL(nday) / REAL( nobis(imois) ) + 0.5 - i15 !!! <== Caution nobis hard coded !!! 324 334 flxnow(:,:,ji) = ( ( 1. - zxy) * flxdta(:,:,ji,1) + zxy * flxdta(:,:,ji,2) ) … … 370 380 IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN 371 381 ! 372 373 374 375 376 382 ! zsst(:,:) = gsst(:,:) / REAL( nfbulk, wp ) * tmask(:,:,1) ! mean sst in K 383 ! zsss(:,:) = gsss(:,:) / REAL( nfbulk ) * tmask(:,:,1) ! mean tn_ice in K 384 385 ! where ( zsst(:,:) .eq. 0.e0 ) zsst(:,:) = rt0 !lb vilain !!??? 386 ! where ( zsss(:,:) .eq. 0.e0 ) zsss(:,:) = rt0 !lb // 377 387 378 388 !!gm better coded: 379 389 ! Caution set to rt0 over land, not 0 Kelvin (otherwise floating point exception in bulk computation) 380 !zcof = 1. / REAL( nfbulk, wp )381 ! zsst(:,:) = gsst(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1) ! mean sst in K382 ! zsss(:,:) = gsss(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1) ! mean tn_ice in K390 zcof = 1. / REAL( nfbulk, wp ) 391 zsst(:,:) = gsst(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1)) ! mean sst in K 392 zsss(:,:) = gsss(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1)) ! mean tn_ice in K 383 393 !!gm 384 394 385 395 zut(:,:) = 0.e0 !!gm not necessary but at least first and last column 386 396 zvt(:,:) = 0.e0 387 397 ! lb 388 398 ! Interpolation of surface current at T-point, zut and zvt : 389 DO ji=2, jpi 390 zut(ji,:) = 0.5*(gu(ji-1,:) + gu(ji,:)) / REAL( nfbulk ) 399 ! DO ji=2, jpi 400 ! zut(ji,:) = 0.5*(gu(ji-1,:) + gu(ji,:)) / REAL( nfbulk ) 401 ! END DO 402 ! ! 403 ! DO jj=2, jpj 404 ! zvt(:,jj) = 0.5*(gv(:,jj-1) + gv(:,jj)) / REAL( nfbulk ) 405 ! END DO 406 !!gm better coded 407 zcof = 0.5 / REAL( nfbulk, wp ) 408 DO jj = 2, jpjm1 409 DO ji = 1, jpi 410 zut(ji,jj) = ( gu(ji-1,jj ) + gu(ji,jj) ) * zcof 411 zvt(ji,jj) = ( gv(ji ,jj-1) + gv(ji,jj) ) * zcof 412 END DO 391 413 END DO 392 !393 DO jj=2, jpj394 zvt(:,jj) = 0.5*(gv(:,jj-1) + gv(:,jj)) / REAL( nfbulk )395 END DO396 !!gm better coded397 ! zcof = 0.5 / REAL( nfbulk, wp )398 ! DO jj = 2, jpjm1399 ! DO ji = fs_2, fs_jpim1 ! vector opt.400 ! zut(ji,jj) = ( gu(ji-1,jj ) + gu(ji,jj) ) * zcof401 ! zvt(ji,jj) = ( gv(ji ,jj-1) + gv(ji,jj) ) * zcof402 ! END DO403 ! END DO404 414 !!gm 405 415 … … 424 434 qlw_ice(:,:) = 0.95 * ( flxnow(:,:,6) - Stef*zsss(:,:)*zsss(:,:)*zsss(:,:)*zsss(:,:) ) ! Long Waves (Infra-red) 425 435 !-------------------------------------------------------------------------------! 426 427 436 428 437 ! ----------------------------------------------------------------------------- ! 429 438 ! II Turbulent FLUXES ! … … 437 446 438 447 !--- Module of relative wind 439 440 448 ! dUnormt(:,:) = sqrt( (flxnow(:,:,2) - zut(:,:))*(flxnow(:,:,2) - zut(:,:)) & 449 ! & + (flxnow(:,:,3) - zvt(:,:))*(flxnow(:,:,3) - zvt(:,:)) ) 441 450 !!gm more efficient coding: 442 ! DO jj = 1, jpi 443 !DO ji = 1, jpi444 !zzu = flxnow(ji,jj,2) - zut(ji,jj)445 !zzv = flxnow(ji,jj,3) - zvt(ji,jj)446 !dUnormt(ji,jj) = SQRT( zzu*zzu + zzv*zzv )447 !END DO448 !END DO451 DO jj = 1, jpj 452 DO ji = 1, jpi 453 zzu = flxnow(ji,jj,2) - zut(ji,jj) 454 zzv = flxnow(ji,jj,3) - zvt(ji,jj) 455 dUnormt(ji,jj) = SQRT( zzu*zzu + zzv*zzv ) 456 END DO 457 END DO 449 458 !!gm 450 ! lb451 !452 459 !--- specific humidity at temp SST 453 460 qsatw(:,:) = 0.98*640380*exp(-5107.4/zsst(:,:))/rhoa … … 462 469 IF( kt == nit000 .AND. lwp ) THEN 463 470 WRITE(numout,*) 464 WRITE(numout,*)' flx_core : global CORE fields in NetCDF format'465 WRITE(numout,*)' ~~~~~~~~~ '466 471 WRITE(numout,*) ' Calling TURB_CORE_2Z for bulk transfert coefficients' 467 WRITE(numout,*)468 472 ENDIF 469 473 !! If air temp. and spec. hum. are given at different height (2m) than wind (10m) : … … 473 477 IF( kt == nit000 .AND. lwp ) THEN 474 478 WRITE(numout,*) 475 WRITE(numout,*)' flx_core : global CORE fields in NetCDF format'476 WRITE(numout,*)' ~~~~~~~~~~ '477 479 WRITE(numout,*) ' Calling TURB_CORE_1Z for bulk transfert coefficients' 478 WRITE(numout,*)479 480 ENDIF 480 481 !! If air temp. and spec. hum. are given at same height than wind (10m) : … … 494 495 CALL lbc_lnk( tauxt(:,:), 'T', -1. ) !!gm seems not decessary at this point.... 495 496 CALL lbc_lnk( tauyt(:,:), 'T', -1. ) 496 ! 497 498 ! Katabatic winds parametrization 499 ! ------------------------------- 500 IF( ln_kata ) THEN 501 ! 502 IF( kt == nit000 ) THEN 503 IF (lwp ) WRITE(numout,*) ' Katabatic winds parametrization is active' 504 clnamekata = 'katamask.nc' 505 clvarkatax = 'katamaskx' 506 clvarkatay = 'katamasky' 507 508 #if defined key_agrif 509 IF ( .NOT. Agrif_Root() ) THEN 510 clnamekata = TRIM(Agrif_CFixed())//'_'//TRIM(clnamekata) 511 ENDIF 512 #endif 513 CALL iom_open( TRIM(clnamekata), inum ) 514 CALL iom_get ( inum, jpdom_data, TRIM(clvarkatax), rmskkatax ) 515 CALL iom_get ( inum, jpdom_data, TRIM(clvarkatay), rmskkatay ) 516 CALL iom_close ( inum ) 517 518 WHERE( (rmskkatax < 0.000001) .AND. (rmskkatax > -0.000001) ) 519 rmskkatax=1 520 rmskkatay=1 521 END WHERE 522 523 CALL lbc_lnk( rmskkatax(:,:), 'T', 1. ) ; CALL lbc_lnk( rmskkatay(:,:), 'T', 1. ) 524 525 IF(MAXVAL(rmskkatax) > 6.00001) THEN 526 WRITE(ctmp1,*) 'Problem in the data mask : maxval = ',MAXVAL(rmskkatax),' ( it is GT 6)' 527 CALL ctl_stop( ctmp1 ) 528 ENDIF 529 ENDIF 530 531 ! apply katabatic wind enhancement on grid T (before projection) 532 tauxt(:,:) = rmskkatax(:,:) * tauxt(:,:) 533 tauyt(:,:) = rmskkatay(:,:) * tauyt(:,:) 534 ! 535 ENDIF 497 536 498 ! 537 499 !Tau_x at U-point … … 555 517 ! --------------------------------- 556 518 ! 557 IF 519 IF( ln_2m ) THEN 558 520 !! 559 521 !! Values of temp. and hum. adjusted to 10m must be used instead of 2m values … … 574 536 END IF 575 537 !! 576 qla_oce(:,:) = Lv * evap(:,:) ! right sign for ocean538 qla_oce(:,:) = -Lv * evap(:,:) ! right sign for ocean 577 539 578 540 ! … … 581 543 ! 582 544 ! Sensible Heat : 583 qsb_ice(:,:) = rhoa*cp*Cice*( zsss(:,:) - flxnow(:,:,7))*dUnormt(:,:) !lb use545 qsb_ice(:,:) = rhoa*cp*Cice*( flxnow(:,:,7) - zsss(:,:) )*dUnormt(:,:) !lb use 584 546 ! !lb dUnormt??? 585 547 ! Latent Heat : !lb or rather Unormt? 586 qla_ice(:,:) = Ls*rhoa*Cice*( qsat(:,:) - flxnow(:,:,4))*dUnormt(:,:)548 qla_ice(:,:) = Ls*rhoa*Cice*( flxnow(:,:,4) - qsat(:,:) )*dUnormt(:,:) 587 549 ! 588 550 !------------------------------------------------------------------------------------- 589 590 591 551 552 592 553 ! ----------------------------------------------------------------------------- ! 593 554 ! III Total FLUXES ! … … 596 557 ! III.1) Downward Non Solar flux over ocean 597 558 ! ----------------------------------------- 598 qnsr_oce(:,:) = qlw_oce(:,:) - qsb_oce(:,:) -qla_oce(:,:)559 qnsr_oce(:,:) = qlw_oce(:,:) + qsb_oce(:,:) + qla_oce(:,:) 599 560 ! 600 561 ! III.1) Downward Non Solar flux over ice 601 562 ! --------------------------------------- 602 qnsr_ice(:,:) = qlw_ice(:,:) - qsb_ice(:,:) -qla_ice(:,:)563 qnsr_ice(:,:) = qlw_ice(:,:) + qsb_ice(:,:) + qla_ice(:,:) 603 564 ! 604 565 !-------------------------------------------------------------------------------! 605 606 607 566 608 567 ! 6. TOTAL NON SOLAR SENSITIVITY 609 568 610 569 dqlw_ice(:,:)= 4.0*0.95*Stef*zsss(:,:)*zsss(:,:)*zsss(:,:) 611 570 612 571 ! d qla_ice/ d zsss 613 572 dqla_ice(:,:) = -Ls*Cice*0.98*11637800/(rhoa*rhoa) & 614 573 * (-5897.8)/(zsss(:,:)*zsss(:,:)) & 615 574 * exp(-5897.8/zsss(:,:)) * dUnormt(:,:) 616 575 617 576 ! d qsb_ice/ d zsss 618 577 dqsb_ice(:,:) = rhoa * cp * Cice * dUnormt(:,:) 619 578 620 579 dqns_ice(:,:) = - ( dqlw_ice(:,:) + dqsb_ice(:,:) + dqla_ice(:,:) ) 621 622 580 623 581 !-------------------------------------------------------------------- 624 582 ! FRACTION of net shortwave radiation which is not absorbed in the 625 583 ! thin surface layer and penetrates inside the ice cover 626 584 ! ( Maykut and Untersteiner, 1971 ; Elbert and Curry, 1993 ) 627 585 628 586 catm1(:,:) = 1.0 - 0.3 ! flxnow(:,:,8) 629 587 fr1_i0(:,:) = & … … 631 589 fr2_i0(:,:) = & 632 590 (0.82 * catm1(:,:) + 0.65 * flxnow(:,:,8)) 633 634 591 592 635 593 ! Total PRECIPITATION (kg/m2/s) 636 tprecip(:,:) = flxnow(:,:,1)594 tprecip(:,:) = alpha_precip*flxnow(:,:,1) 637 595 ! rename precipitation for freshwater budget calculations: 638 watm(:,:) = flxnow(:,:,1)*1000 639 ! 640 596 ! watm(:,:) = flxnow(:,:,1)*1000 597 watm(:,:) = flxnow(:,:,1)*rday 598 ! 599 641 600 ! SNOW PRECIPITATION (kg/m2/s) 642 sprecip(:,:) = flxnow(:,:,8)643 601 sprecip(:,:) = alpha_precip*flxnow(:,:,8) 602 644 603 !--------------------------------------------------------------------- 645 646 647 604 648 605 CALL lbc_lnk( taux (:,:) , 'U', -1. ) 649 606 CALL lbc_lnk( tauy (:,:) , 'V', -1. ) … … 679 636 gu (:,:) = 0.e0 680 637 gv (:,:) = 0.e0 681 638 682 639 ! Degugging print (A.M. Treguier) 683 640 ! write (55) taux, tauy, qnsr_oce, qsb_ice, qnsr_ice, qla_ice, dqns_ice & 684 641 ! & , dqla_ice, fr1_i0, fr2_i0, qlw_oce, qla_oce, qsb_oce, evap 685 642 ! write(numout,*) 'written 14 arrays on unit fort.55' 686 687 643 644 688 645 ENDIF 689 646 690 647 ! ------------------- ! 691 648 ! Last call kt=nitend ! 692 649 ! ------------------- ! 693 650 694 651 ! Closing of the numflx file (required in mpp) 695 652 696 653 IF( kt == nitend ) THEN 697 654 DO ji=1, jpfile … … 699 656 END DO 700 657 ENDIF 701 658 IF ( lrst_oce ) CALL core_rst( kt, 'WRITE' ) 659 702 660 END SUBROUTINE flx 703 704 705 SUBROUTINE TURB_CORE_1Z( zzu, T_0, T_a, q_sat, q_a, & 706 & dU , C_d, C_h , C_e ) 661 662 SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a, & 663 & dU, Cd, Ch, Ce ) 707 664 !!---------------------------------------------------------------------- 708 665 !! *** ROUTINE turb_core *** 709 666 !! 710 !! ** Purpose : Computes turbulent transfert coefficients of surface 711 !! fluxes according to Large & Yeager (2004) .667 !! ** Purpose : Computes turbulent transfert coefficients of surface 668 !! fluxes according to Large & Yeager (2004) 712 669 !! 713 670 !! ** Method : I N E R T I A L D I S S I P A T I O N M E T H O D … … 721 678 !! History : 722 679 !! ! XX-XX (??? ) Original code 723 !! 9.0 ! 05-08 (L. Brodeau) Optimisation680 !! 9.0 ! 05-08 (L. Brodeau) Rewriting and optimization 724 681 !!---------------------------------------------------------------------- 725 682 !! * Arguments 726 REAL(wp), INTENT( in ) :: & 727 zzu ! height for all given atmospheric variables[m]728 REAL(wp), INTENT( in), DIMENSION(jpi,jpj) :: &729 T_0, & ! sea surface temperature[Kelvin]730 T_a, & ! potential air temperature at zzu[Kelvin]731 q_sat, & ! sea surface specific humidity at zzu[kg/kg]732 q_a, & ! specific air humidity at zzu[kg/kg]733 dU ! relative wind module |U(zzu)-U(0)| at zzu[m/s]734 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: &735 C _d, & ! transfercoefficient for momentum (tau)736 C _h, & ! transfer coefficient for sensible heat(Q_sens)737 C _e ! tansfert coefficient for evaporation(Q_lat)683 684 REAL(wp), INTENT(in) :: zu ! altitude of wind measurement [m] 685 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: & 686 sst, & ! sea surface temperature [Kelvin] 687 T_a, & ! potential air temperature [Kelvin] 688 q_sat, & ! sea surface specific humidity [kg/kg] 689 q_a, & ! specific air humidity [kg/kg] 690 dU ! wind module |U(zu)-U(0)| [m/s] 691 REAL(wp), intent(out), DIMENSION(jpi,jpj) :: & 692 Cd, & ! transfert coefficient for momentum (tau) 693 Ch, & ! transfert coefficient for temperature (Q_sens) 694 Ce ! transfert coefficient for evaporation (Q_lat) 738 695 739 696 !! * Local declarations … … 745 702 Ce_n10, & ! 10m neutral latent coefficient 746 703 Ch_n10, & ! 10m neutral sensible coefficient 747 Cd, & ! drag coefficient748 Ce, & ! latent coefficient749 Ch, & ! sensible coefficient750 704 sqrt_Cd_n10, & ! root square of Cd_n10 751 705 sqrt_Cd, & ! root square of Cd … … 755 709 U_star, & ! turb. scale of velocity fluct. 756 710 L, & ! Monin-Obukov length [m] 757 zeta, & ! stability parameter at height zzu 758 X2, X, & 759 psi_m, & 760 psi_h, & 761 U_n10, & ! neutral wind velocity at 10m [m] 762 xlogt 763 764 INTEGER :: jk ! doomy loop for iterations 765 INTEGER, PARAMETER :: nit = 3 ! number of iterations 766 711 zeta, & ! stability parameter at height zu 712 U_n10, & ! neutral wind velocity at 10m [m] 713 xlogt, xct, zpsi_h, zpsi_m 714 !! 715 INTEGER :: j_itt 716 INTEGER, PARAMETER :: nb_itt = 3 767 717 INTEGER, DIMENSION(jpi,jpj) :: & 768 stab, & ! stability test integer 769 stabit ! stability within iterative loop 770 771 REAL(wp), PARAMETER :: & 772 pi = 3.14159, & ! Pi 773 grav = 9.8, & ! gravity 774 kappa = 0.4 ! von Karman's constant 775 !!---------------------------------------------------------------------- 776 !! 777 !! ------------- 778 !! S T A R T 779 !! ------------- 780 !! 781 !! I.1 ) Preliminary stuffs 782 !! ------------------------ 783 !! 718 stab ! 1st guess stability test integer 719 720 REAL(wp), PARAMETER :: & 721 grav = 9.8, & ! gravity 722 kappa = 0.4 ! von Karman s constant 723 724 !! * Start 784 725 !! Air/sea differences 785 !! -------------------786 726 dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s 787 dT = T_a - T_0! assuming that T_a is allready the potential temp. at zzu788 dq 727 dT = T_a - sst ! assuming that T_a is allready the potential temp. at zzu 728 dq = q_a - q_sat 789 729 !! 790 730 !! Virtual potential temperature 791 !! -----------------------------792 731 T_vpot = T_a*(1. + 0.608*q_a) 793 732 !! 794 !! 795 !! I.2 ) Computing Neutral Drag Coefficient 796 !! ---------------------------------------- 797 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! \\ L & Y eq. (6a) 733 !! Neutral Drag Coefficient 734 stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 735 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 798 736 sqrt_Cd_n10 = sqrt(Cd_n10) 799 !! 800 Ce_n10 = 1E-3 * ( 34.6 * sqrt_Cd_n10 ) ! \\ L & Y eq. (6b) 801 !! 802 !! First guess of stabilitty : 803 stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 804 Ch_n10 = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) ! \\ L & Y eq. (6c), (6d) 805 !! 806 !! Initializing transfert coefficients with their first guess neutral equivalents : 807 Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 808 !! 809 !! 810 !! 811 !! II ) Now starting iteration loop (IDM) 812 !! ---------------------------------------- 813 !! 814 DO jk=1, nit 815 !! 816 sqrt_Cd = sqrt(Cd) 817 !! 818 !! Turbulent scales : 819 !! ------------------ 820 U_star = sqrt_Cd*dU10 ! \\ L & Y eq. (7a) 821 T_star = Ch/sqrt_Cd*dT ! \\ L & Y eq. (7b) 822 q_star = Ce/sqrt_Cd*dq ! \\ L & Y eq. (7c) 823 !! 824 !! Estimate the Monin-Obukov length : 825 !! ---------------------------------- 826 !dbg1 = T_star/T_vpot 827 !dbg2 = q_star/(q_a + 1./0.608) 828 L = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) ) 829 !! 830 !! Stability parameters : 831 !! ---------------------- 832 zeta = zzu/L 833 zeta = sign( min(abs(zeta),10.0), zeta ) 834 !! 835 !! Psis, L & Y eq. (8c), (8d), (8e) : 836 !! ---------------------------------- 837 X2 = sqrt(abs(1. - 16.*zeta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) 838 !! 839 stabit = 0.5 + sign(0.5,zeta) 840 !! 841 ! dbg1 = -5*zeta*stabit 842 ! dbg2 = 2*log((1. + X)/2) 843 ! dbg3 = log((1. + X2)/2) 844 ! dbg4 = 2*atan(X) 845 846 psi_m = -5*zeta*stabit & ! Stable 847 & + (1 - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 848 !! 849 psi_h = -5*zeta*stabit & ! Stable 850 & + (1 - stabit)*(2*log( (1. + X2)/2 )) ! Unstable 851 !! 852 !! 853 !! Shifting the wind speed to 10m and neutral stability : 854 !! ------------------------------------------------------ 855 U_n10 = dU10*1./(1. + sqrt(Cd_n10)/kappa*(log(zzu/10.) - psi_m)) ! \\ L & Y eq. (9a) 856 !! 857 !! 858 !! Updating the neutral 10m transfer coefficients : 859 !! ------------------------------------------------ 860 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! \\ L & Y eq. (6a) 861 sqrt_Cd_n10 = sqrt(Cd_n10) 862 !! 863 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! \\ L & Y eq. (6b) 864 !! 865 stab = 0.5 + sign(0.5,zeta) 866 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! \\ L & Y eq. (6c), (6d) 867 !! 868 !! 869 !! Shifting the neutral 10m transfer coefficients to ( zzu , zeta ) : 870 !! -------------------------------------------------------------------- 871 !! Problem here, formulation used within L & Y differs from the one provided 872 !! in their fortran code (only for Ce and Ch) 873 !! 874 Cd = Cd_n10/(1. + sqrt_Cd_n10/kappa*(log(zzu/10) - psi_m))**2 ! \\ L & Y eq. (10a) 875 !! 876 xlogt = log(zzu/10) - psi_h 877 !? Ch = Ch_n10*sqrt(Cd/Cd_n10)/(1. + Ch_n10/(kappa*sqrt_Cd_n10)*xlogt) 878 Ch = Ch_n10/( 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 )**2 ! \\ L & Y eq. (10b) 879 !! 880 !? Ce = Ce_n10*sqrt(Cd/Cd_n10)/(1. + Ce_n10/(kappa*sqrt_Cd_n10)*xlogt) 881 Ce = Ce_n10/( 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 )**2 ! \\ L & Y eq. (10c) 882 !! 883 !! 884 END DO 885 !! 886 C_d = Cd 887 C_h = Ch 888 C_e = Ce 889 !! 890 END SUBROUTINE TURB_CORE_1Z 891 892 893 SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu) 737 Ce_n10 = 1E-3 * ( 34.6 * sqrt_Cd_n10 ) ! L & Y eq. (6b) 738 Ch_n10 = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) ! L & Y eq. (6c), (6d) 739 !! 740 !! Initializing transfert coefficients with their first guess neutral equivalents : 741 Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd) 742 743 !! * Now starting iteration loop 744 DO j_itt=1, nb_itt 745 !! Turbulent scales : 746 U_star = sqrt_Cd*dU10 ! L & Y eq. (7a) 747 T_star = Ch/sqrt_Cd*dT ! L & Y eq. (7b) 748 q_star = Ce/sqrt_Cd*dq ! L & Y eq. (7c) 749 750 !! Estimate the Monin-Obukov length : 751 L = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) ) 752 753 !! Stability parameters : 754 zeta = zu/L ; zeta = sign( min(abs(zeta),10.0), zeta ) 755 zpsi_h = psi_h(zeta) 756 zpsi_m = psi_m(zeta) 757 758 !! Shifting the wind speed to 10m and neutral stability : 759 U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) ! L & Y eq. (9a) 760 761 !! Updating the neutral 10m transfer coefficients : 762 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 763 sqrt_Cd_n10 = sqrt(Cd_n10) 764 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 765 stab = 0.5 + sign(0.5,zeta) 766 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c), (6d) 767 768 !! Shifting the neutral 10m transfer coefficients to ( zu , zeta ) : 769 !! 770 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m) 771 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 772 !! 773 xlogt = log(zu/10.) - zpsi_h 774 !! 775 xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 776 Ch = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct 777 !! 778 xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 779 Ce = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 780 !! 781 END DO 782 !! 783 END SUBROUTINE TURB_CORE_1Z 784 785 786 SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu) 894 787 !!---------------------------------------------------------------------- 895 788 !! *** ROUTINE turb_core *** … … 907 800 !! Large & Yeager, 2004 : ??? 908 801 !! History : 909 !! ! XX-XX (??? ) Original code 910 !! 9.0 ! 05-08 (L. Brodeau) Optimisation 802 !! 9.0 ! 06-12 (L. Brodeau) Original code for 2Z 911 803 !!---------------------------------------------------------------------- 912 804 !! * Arguments 913 REAL(wp), INTENT( in ) :: & 914 zt, & ! height for T_zt and q_zt [m] 915 zu ! height for dU [m] 916 !! 805 REAL(wp), INTENT(in) :: & 806 zt, & ! height for T_zt and q_zt [m] 807 zu ! height for dU [m] 917 808 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: & 918 sst, & ! sea surface temperature [Kelvin] 919 T_zt, & ! potential air temperature [Kelvin] 920 q_sat, & ! sea surface specific humidity [kg/kg] 921 q_zt, & ! specific air humidity [kg/kg] 922 dU ! wind module |U(zu)-U(0)| [m/s] 923 !! 809 sst, & ! sea surface temperature [Kelvin] 810 T_zt, & ! potential air temperature [Kelvin] 811 q_sat, & ! sea surface specific humidity [kg/kg] 812 q_zt, & ! specific air humidity [kg/kg] 813 dU ! relative wind module |U(zu)-U(0)| [m/s] 924 814 REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: & 925 Cd, Ch, Ce, & 926 T_zu, & ! air temp. shifted at zu [K] 927 q_zu ! spec. hum. shifted at zu [kg/kg] 815 Cd, & ! transfer coefficient for momentum (tau) 816 Ch, & ! transfer coefficient for sensible heat (Q_sens) 817 Ce, & ! transfert coefficient for evaporation (Q_lat) 818 T_zu, & ! air temp. shifted at zu [K] 819 q_zu ! spec. hum. shifted at zu [kg/kg] 928 820 929 821 !! * Local declarations 930 REAL(wp), DIMENSION(jpi,jpj) :: & 931 dU10, & ! dU [m/s] 932 dT, & ! air/sea temperature differeence [K] 933 dq, & ! air/sea humidity difference [K] 934 Cd_n10, & ! 10m neutral drag coefficient 935 Ce_n10, & ! 10m neutral latent coefficient 936 Ch_n10, & ! 10m neutral sensible coefficient 937 sqrt_Cd_n10, & ! root square of Cd_n10 938 sqrt_Cd, & ! root square of Cd 939 T_vpot_u, & ! virtual potential temperature [K] 940 T_star, & ! turbulent scale of tem. fluct. 941 q_star, & ! turbulent humidity of temp. fluct. 942 U_star, & ! turb. scale of velocity fluct. 943 L, & ! Monin-Obukov length [m] 944 zeta_u, & ! stability parameter at height zu 945 zeta_t, & ! stability parameter at height zt 946 U_n10, & ! neutral wind velocity at 10m [m] 947 xlogt, xct 948 !! 822 REAL(wp), DIMENSION(jpi,jpj) :: & 823 dU10, & ! dU [m/s] 824 dT, & ! air/sea temperature differeence [K] 825 dq, & ! air/sea humidity difference [K] 826 Cd_n10, & ! 10m neutral drag coefficient 827 Ce_n10, & ! 10m neutral latent coefficient 828 Ch_n10, & ! 10m neutral sensible coefficient 829 sqrt_Cd_n10, & ! root square of Cd_n10 830 sqrt_Cd, & ! root square of Cd 831 T_vpot_u, & ! virtual potential temperature [K] 832 T_star, & ! turbulent scale of tem. fluct. 833 q_star, & ! turbulent humidity of temp. fluct. 834 U_star, & ! turb. scale of velocity fluct. 835 L, & ! Monin-Obukov length [m] 836 zeta_u, & ! stability parameter at height zu 837 zeta_t, & ! stability parameter at height zt 838 U_n10, & ! neutral wind velocity at 10m [m] 839 xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m 840 841 INTEGER :: j_itt 949 842 INTEGER, PARAMETER :: nb_itt = 3 ! number of itterations 950 !!951 !! Some physical constants :843 INTEGER, DIMENSION(jpi,jpj) :: & 844 & stab ! 1st stability test integer 952 845 REAL(wp), PARAMETER :: & 953 grav = 9.8, & ! gravity 954 kappa = 0.4 ! von Karman's constant 955 !! 956 INTEGER :: j_itt 957 INTEGER, DIMENSION(jpi,jpj) :: & 958 stab ! 1st stability test integer 959 !!---------------------------------------------------------------------- 960 !! 961 !! ------------- 962 !! S T A R T 963 !! ------------- 964 !! 965 !! I.1 ) Preliminary stuffs 966 !! ------------------------ 967 !! 846 grav = 9.8, & ! gravity 847 kappa = 0.4 ! von Karman's constant 848 849 !! * Start 850 968 851 !! Initial air/sea differences 969 dU10 = max(0.5, dU) ; dT = T_zt - sst ; dq = q_zt - q_sat 970 !! 852 dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s 853 dT = T_zt - sst 854 dq = q_zt - q_sat 855 971 856 !! Neutral Drag Coefficient : 972 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE857 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE 973 858 Cd_n10 = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 974 859 sqrt_Cd_n10 = sqrt(Cd_n10) 975 860 Ce_n10 = 1E-3*( 34.6 * sqrt_Cd_n10 ) 976 861 Ch_n10 = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab)) 977 !! 978 !! I.2 ) Computing Neutral Drag Coefficient 979 !! ---------------------------------------- 862 980 863 !! Initializing transf. coeff. with their first guess neutral equivalents : 981 864 Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd) 982 !! 865 983 866 !! Initializing z_u values with z_t values : 984 867 T_zu = T_zt ; q_zu = q_zt 985 !! 986 !! II ) Now starting iteration loop (IDM) 987 !! ---------------------------------------- 988 !! 868 869 !! * Now starting iteration loop 989 870 DO j_itt=1, nb_itt 990 !! 991 !! Updating air/sea differences : 992 dT = T_zu - sst ; dq = q_zu - q_sat 993 !! 994 !! Updating virtual potential temperature at zu : 995 T_vpot_u = T_zu*(1. + 0.608*q_zu) 996 !! 997 !! Updating turbulent scales : (L & Y eq. (7)) 998 U_star = sqrt_Cd*dU10 ; T_star = Ch/sqrt_Cd*dT ; q_star = Ce/sqrt_Cd*dq 999 !! 1000 !! Estimate the Monin-Obukov length at height zu : 1001 L = (U_star*U_star) & 871 dT = T_zu - sst ; dq = q_zu - q_sat ! Updating air/sea differences 872 T_vpot_u = T_zu*(1. + 0.608*q_zu) ! Updating virtual potential temperature at zu 873 U_star = sqrt_Cd*dU10 ! Updating turbulent scales : (L & Y eq. (7)) 874 T_star = Ch/sqrt_Cd*dT ! 875 q_star = Ce/sqrt_Cd*dq ! 876 !! 877 L = (U_star*U_star) & ! Estimate the Monin-Obukov length at height zu 1002 878 & / (kappa*grav/T_vpot_u*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star)) 1003 !!1004 879 !! Stability parameters : 1005 880 zeta_u = zu/L ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 1006 881 zeta_t = zt/L ; zeta_t = sign( min(abs(zeta_t),10.0), zeta_t ) 882 zpsi_hu = psi_h(zeta_u) 883 zpsi_ht = psi_h(zeta_t) 884 zpsi_m = psi_m(zeta_u) 1007 885 !! 1008 886 !! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a)) 1009 U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))) 887 ! U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))) 888 U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) 1010 889 !! 1011 890 !! Shifting temperature and humidity at zu : (L & Y eq. (9b-9c)) 1012 T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 1013 q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 891 ! T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 892 T_zu = T_zt - T_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht) 893 ! q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 894 q_zu = q_zt - q_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht) 1014 895 !! 1015 896 !! q_zu cannot have a negative value : forcing 0 … … 1025 906 !! 1026 907 !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) : 1027 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)) 908 ! xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)) 909 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m) 1028 910 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 1029 911 !! 1030 xlogt = log(zu/10.) - psi_h(zeta_u) 912 ! xlogt = log(zu/10.) - psi_h(zeta_u) 913 xlogt = log(zu/10.) - zpsi_hu 1031 914 !! 1032 915 xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 … … 1039 922 END DO 1040 923 !! 1041 END SUBROUTINE TURB_CORE_2Z 1042 1043 FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 1044 REAL(wp), PARAMETER :: pi = 3.14159 1045 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 1046 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 1047 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 1048 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) 1049 stabit = 0.5 + sign(0.5,zta) 1050 psi_m = -5.*zta*stabit & ! Stable 1051 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 1052 END FUNCTION psi_m 1053 1054 FUNCTION psi_h(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 1055 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 1056 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 1057 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 1058 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2) 1059 stabit = 0.5 + sign(0.5,zta) 1060 psi_h = -5.*zta*stabit & ! Stable 1061 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 1062 END FUNCTION psi_h 924 END SUBROUTINE TURB_CORE_2Z 925 926 FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 927 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 928 929 REAL(wp), PARAMETER :: pi = 3.14159 930 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 931 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 932 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) 933 stabit = 0.5 + sign(0.5,zta) 934 psi_m = -5.*zta*stabit & ! Stable 935 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 936 END FUNCTION psi_m 937 938 FUNCTION psi_h(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 939 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 940 941 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 942 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 943 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2) 944 stabit = 0.5 + sign(0.5,zta) 945 psi_h = -5.*zta*stabit & ! Stable 946 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 947 END FUNCTION psi_h 1063 948 1064 949 … … 1083 968 !! * Modules used 1084 969 USE ice ! ??? 1085 970 1086 971 !! * Arguments 1087 972 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: & … … 1090 975 palbp , & ! albedo of ice under clear sky 1091 976 palcnp ! albedo of ocean under clear sky 1092 977 1093 978 !! * Local variables 1094 979 INTEGER :: & … … 1106 991 zzero = 0.0 , & 1107 992 zone = 1.0 1108 993 1109 994 REAL(wp) :: & 1110 995 zmue14 , & ! zmue**1.4 … … 1129 1014 1130 1015 !!--------------------------------------------------------------------- 1131 1016 1132 1017 !------------------------- 1133 1018 ! Computation of zficeth 1134 1019 !-------------------------- 1135 1020 1136 1021 llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice ) 1137 1022 WHERE ( llmask ) ! ice free of snow and melts … … 1140 1025 zalbfz = alphdi 1141 1026 END WHERE 1142 1027 1143 1028 DO jj = 1, jpj 1144 1029 DO ji = 1, jpi … … 1156 1041 END DO 1157 1042 END DO 1158 1043 1159 1044 !----------------------------------------------- 1160 1045 ! Computation of the snow/ice albedo system 1161 1046 !-------------------------- --------------------- 1162 1047 1163 1048 ! Albedo of snow-ice for clear sky. 1164 1049 !----------------------------------------------- … … 1176 1061 ( albice + hsnif(ji,jj) * ( alphc - albice ) / c2 ) & 1177 1062 & + zihsc2 * alphc 1178 1063 1179 1064 zitmlsn = MAX ( zzero , SIGN ( zone , sist(ji,jj) - rt0_snow ) ) 1180 1065 zalbpsn = zitmlsn * zalbpsnf + ( 1.0 - zitmlsn ) * zalbpsnm 1181 1066 1182 1067 ! Case of ice free of snow. 1183 1068 zalbpic = zficeth(ji,jj) 1184 1069 1185 1070 ! albedo of the system 1186 1071 zithsn = 1.0 - MAX ( zzero , SIGN ( zone , - hsnif(ji,jj) ) ) … … 1188 1073 END DO 1189 1074 END DO 1190 1075 1191 1076 ! Albedo of snow-ice for overcast sky. 1192 1077 !---------------------------------------------- 1193 1078 palb(:,:) = palbp(:,:) + cgren 1194 1079 1195 1080 !-------------------------------------------- 1196 1081 ! Computation of the albedo of the ocean 1197 1082 !-------------------------- ----------------- 1198 1199 1083 1084 1200 1085 ! Parameterization of Briegled and Ramanathan, 1982 1201 1086 zmue14 = zmue**1.4 1202 1087 palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 ) 1203 1088 1204 1089 ! Parameterization of Kondratyev, 1969 and Payne, 1972 1205 1090 palcn(:,:) = 0.06 1206 1091 1207 1092 END SUBROUTINE flx_blk_albedo 1208 1093 1209 1094 SUBROUTINE core_rst( kt, cdrw ) 1095 !!--------------------------------------------------------------------- 1096 !! *** ROUTINE ts_rst *** 1097 !! 1098 !! ** Purpose : Read or write CORE specific variables ( gsss, gu, gv ) 1099 !! 1100 !! ** Method : 1101 !! 1102 !!---------------------------------------------------------------------- 1103 USE iom ! move to top of flxmod to avoid duplication 1104 USE restart ! move to top of flxmod to avoid duplication 1105 ! 1106 INTEGER , INTENT(in) :: kt ! ocean time-step 1107 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1108 ! 1109 !!---------------------------------------------------------------------- 1110 ! 1111 IF( TRIM(cdrw) == 'READ' ) THEN 1112 IF( ln_rstart ) THEN 1113 CALL iom_get( numror, jpdom_local, 'gsss', gsss ) 1114 CALL iom_get( numror, jpdom_local, 'gu', gu ) 1115 CALL iom_get( numror, jpdom_local, 'gv', gv ) 1116 ENDIF ! gsss, gu, gv are initialized in the routine ( may be changed ) 1117 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 1118 CALL iom_rstput( kt, nitrst, numrow, 'gsss', gsss ) 1119 CALL iom_rstput( kt, nitrst, numrow, 'gu', gu ) 1120 CALL iom_rstput( kt, nitrst, numrow, 'gv', gv ) 1121 ENDIF 1122 ! 1123 END SUBROUTINE core_rst
Note: See TracChangeset
for help on using the changeset viewer.
