Index: /dynamico_lmdz/simple_physics/install_param.sh =================================================================== --- /dynamico_lmdz/simple_physics/install_param.sh (revision 4176) +++ /dynamico_lmdz/simple_physics/install_param.sh (revision 4176) @@ -0,0 +1,68 @@ +#!/bin/bash + + +function cmd_install_lmdz() +{ +# Installation du modele en mode de parallelisation mixte MIPxOpenMP + echo "cmd_install_lmdz" + NETCDF_ROOT=$(nc-config --prefix) + echo "NETCDF root : $NETCDF_ROOT" + wget http://www.lmd.jussieu.fr/~lmdz/pub/install_lmdz.sh + export LANG=C # fixes issue with sed on MaxOSX + sed -e 's/veget=1/veget=0/g' install_lmdz.sh > install_lmdz_patched.sh + sed -i -e bak 's/makelmdz_fcm/echo makelmdz/g' install_lmdz_patched.sh + chmod +x install_lmdz_patched.sh + ./install_lmdz_patched.sh -parallel none -v $version +} + +function cmd_get_phyparam() +{ + echo "cmd_get_phyparam" + wget http://www.lmd.jussieu.fr/~hourdin/PARAM/phyparam.20191009.tar + tar xvf phyparam.20191009.tar +} + +function cmd_patch_lmdz() +{ +# Modification du code source pour prendre en compte la physique +# a 20 parametres + echo "cmd_patch_lmdz" + cd $LMDZ/libf + rm -rf phyparam dynphy_lonlat/phyparam + mkdir phyparam dynphy_lonlat/phyparam + cd phyparam + ln -s ../phydev/* . + ln -sf $ROOT/phyparam/* . + ln -sf $ROOT/phyparam/param/* . + cd ../dynphy_lonlat/phyparam + ln -s ../phydev/* . + ln -sf $ROOT/dynphy_lonlat/phyparam/* . + cd $LMDZ + ./makelmdz_fcm -rrtm false -v false -arch local -j 8 -p param -d 32x32x39 gcm +} + + +function cmd_get_param() +{ +# Recupération des parametres d'entree + echo "cmd_get_param" + cd $LMDZ + wget http://www.lmd.jussieu.fr/~hourdin/PARAM/test_param.tar + tar xvf test_param.tar +} + +function cmd_() +{ +# cmd_get_phyparam + cmd_install_lmdz + cmd_patch_lmdz + cmd_get_param + echo puis lancer gcm.e sur TEST_PARAM +} + +# On peut choisir la version de LMDZ a insitaller +version=20191106.trunk +ROOT=$(pwd) +LMDZ=$ROOT/LMDZ$version/modipsl/modeles/LMDZ + +cmd_$1 Index: /dynamico_lmdz/simple_physics/phyparam/param/callkeys.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/callkeys.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/callkeys.h (revision 4176) @@ -0,0 +1,6 @@ + COMMON/callkeys/callrad,calldifv,calladj,callcond,callsoil, + s season,diurnal,lverbose,iradia,period_sort + LOGICAL callrad,calldifv,calladj,callcond,callsoil, + s season,diurnal,lverbose + INTEGER iradia + real period_sort Index: /dynamico_lmdz/simple_physics/phyparam/param/clesphys.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/clesphys.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/clesphys.h (revision 4176) @@ -0,0 +1,44 @@ +! +! $Header: /home/cvsroot/LMDZ4/libf/phylmd/clesphys.h,v 1.5 2006/05/09 09:46:14 fairhead Exp $ +! +c..include cles_phys.h +c + LOGICAL cycle_diurne,soil_model,new_oliq,ok_orodr,ok_orolf + LOGICAL ok_limitvrai + INTEGER nbapp_rad, iflag_con + REAL co2_ppm, solaire + REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12 + REAL*8 CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt +cIM simulateur ISCCP + INTEGER top_height, overlap +cIM seuils cdrm, cdrh + REAL cdmmax, cdhmax +cIM param. stabilite s/ terres et en dehors + REAL ksta, ksta_ter +cIM ok_kzmin : clef calcul Kzmin dans la CL de surface cf FH + LOGICAL ok_kzmin +cIM lev_histhf : niveau sorties 6h +cIM lev_histday : niveau sorties journalieres +cIM lev_histmth : niveau sorties mensuelles + INTEGER lev_histhf, lev_histday, lev_histmth + CHARACTER*4 type_run + LOGICAL ok_isccp, ok_regdyn + REAL lonmin_ins, lonmax_ins, latmin_ins, latmax_ins + REAL ecrit_ins, ecrit_hf, ecrit_hf2mth, ecrit_day + REAL ecrit_mth, ecrit_tra, ecrit_reg + REAL freqin_isccp, freqout_isccp + INTEGER :: ip_ebil_phy + LOGICAL ok_slab_sicOBS + + COMMON/clesphys/cycle_diurne, soil_model, new_oliq, + S ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad, iflag_con + S , co2_ppm, solaire, RCO2, RCH4, RN2O, RCFC11, RCFC12 + S , CH4_ppb, N2O_ppb, CFC11_ppt, CFC12_ppt + S , top_height, overlap, cdmmax, cdhmax, ksta, ksta_ter + S , ok_kzmin, lev_histhf, lev_histday, lev_histmth + S , type_run, ok_isccp, ok_regdyn + S , lonmin_ins, lonmax_ins, latmin_ins, latmax_ins + S , ecrit_ins, ecrit_hf, ecrit_hf2mth, ecrit_day + S , ecrit_mth, ecrit_tra, ecrit_reg + S , freqin_isccp, freqout_isccp, ip_ebil_phy + S , ok_slab_sicOBS Index: /dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/coefdifv.F (revision 4176) @@ -0,0 +1,78 @@ + SUBROUTINE coefdifv( np,nlev, + $ pu,pv,ph,pphi, + $ pcdzv,pcdzh) + USE constlim + IMPLICIT NONE + +c======================================================================= +c +c calcul des coeficients de la diffusion verticale +c +c======================================================================= +c +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "dimensions.h" +#include "comcstfi.h" + +c Arguments: +c ---------- + + INTEGER np,nlev + REAL pu(np,nlev),pv(np,nlev) + REAL ph(np,nlev),pphi(np,nlev) + REAL pcdzv(np,nlev),pcdzh(np,nlev) + +c Local: +c ------ + + INTEGER ilev,ip + REAL z1(np),z2(np) + REAL z1odz(np),zdzu2(np),zhbar(np) + +c----------------------------------------------------------------------- +c initialisations: +c ---------------- + +c----------------------------------------------------------------------- +c couche de surface: +c ------------------ + + DO 210 ip=1,np + z1(ip)=pu(ip,1)*pu(ip,1)+pv(ip,1)*pv(ip,1) + pcdzv(ip,1)=dgcdrag(ip)*(SQRT(z1(ip))+1.) / ph(ip,1) + pcdzh(ip,1)=pcdzv(ip,1) +210 CONTINUE +c +c +c----------------------------------------------------------------------- +c autres couches: +c --------------- + + PRINT*,'Coeff k' + DO 310 ilev=1,nlev-1 + DO 320 ip=1,np + z1(ip)=pu(ip,ilev+1)-pu(ip,ilev) + z2(ip)=pv(ip,ilev+1)-pv(ip,ilev) + z1(ip)=z1(ip)*z1(ip)+z2(ip)*z2(ip) + z1odz(ip)=g/(pphi(ip,ilev+1)-pphi(ip,ilev)) + zdzu2(ip)=z1(ip)*z1odz(ip)*z1odz(ip) + zhbar(ip)=0.5*(ph(ip,ilev)+ph(ip,ilev+1)) + z1(ip)=( (ph(ip,ilev+1)-ph(ip,ilev))*z1odz(ip)- + $ cpgam ) * ccdzh/zhbar(ip) + pcdzv(ip,ilev+1)=cdzconst(ilev+1)* + $ SQRT(AMAX1(zdzu2(ip)-z1(ip),cdzmin)) / +c $ SQRT(cdzmin) / + $ (zhbar(ip)*zhbar(ip)) + IF(ip.eq.np/2+1) PRINT*,lmixmin*lmixmin* + s SQRT(AMAX1(zdzu2(ip)-z1(ip),cdzmin)) + pcdzh(ip,ilev+1)=pcdzv(ip,ilev+1) +320 CONTINUE +310 CONTINUE + +c----------------------------------------------------------------------- + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/comcstfi.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/comcstfi.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/comcstfi.h (revision 4176) @@ -0,0 +1,9 @@ +c----------------------------------------------------------------------- +c INCLUDE comcstfi.h + + COMMON/comcstfi/ + * pi,rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg + + REAL pi,rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg + +c----------------------------------------------------------------------- Index: /dynamico_lmdz/simple_physics/phyparam/param/comcstphy.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/comcstphy.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/comcstphy.h (revision 4176) @@ -0,0 +1,2 @@ + real rradius,rr,rg,rcpp + common/comcstphy/rradius,rr,rg,rcpp Index: /dynamico_lmdz/simple_physics/phyparam/param/comgeomfi.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/comgeomfi.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/comgeomfi.F90 (revision 4176) @@ -0,0 +1,41 @@ +module comgeomfi + + real,save,allocatable :: long(:) + real,save,allocatable :: lati(:) + real,save,allocatable :: area(:) + real,save,allocatable :: sinlon(:) + real,save,allocatable :: coslon(:) + real,save,allocatable :: sinlat(:) + real,save,allocatable :: coslat(:) + real,save :: totarea + integer,save :: ngridmax,nlayermx,nsoilmx +!$OMP THREADPRIVATE(long,lati,area,sinlon,coslon,sinlat,coslat,totarea) +!$OMP THREADPRIVATE(ngridmax,nlayermx,nsoilmx) +contains + + subroutine InitComgeomfi + USE mod_phys_lmdz_para + USE dimphy, ONLY : klon,klev + USE geometry_mod, ONLY : latitude_deg,longitude_deg + implicit none + + + print*,'Dans initcomgeomfi ' + ngridmax=klon_omp + nlayermx=klev + nsoilmx=10 + print*,'ngridmax,nlayermx',ngridmax,nlayermx + + allocate(long(klon_omp)) + allocate(lati(klon_omp)) + long=longitude_deg + lati=latitude_deg + allocate(area(klon_omp)) + allocate(sinlon(klon_omp)) + allocate(coslon(klon_omp)) + allocate(sinlat(klon_omp)) + allocate(coslat(klon_omp)) + + end subroutine InitComgeomfi + +end module comgeomfi Index: /dynamico_lmdz/simple_physics/phyparam/param/comlw.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/comlw.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/comlw.h (revision 4176) @@ -0,0 +1,8 @@ +c----------------------------------------------------------------------- +c INCLUDE 'comlw.h' + + REAL teq(nlayermx),dtlw(nlayermx),taulw(nlayermx) + + COMMON/comlw/teq,dtlw,taulw + +c----------------------------------------------------------------------- Index: /dynamico_lmdz/simple_physics/phyparam/param/comorbit.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/comorbit.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/comorbit.h (revision 4176) @@ -0,0 +1,12 @@ +c----------------------------------------------------------------------- +c INCLUDE 'comorbit.h' + + COMMON/comorbit/aphelie,periheli,year_day, + $ peri_day,timeperi,obliquit, + $ e_elips,p_elips,unitastr,pi + + REAL aphelie,periheli,year_day, + $ peri_day,timeperi,obliquit, + $ e_elips,p_elips,unitastr,pi + +c----------------------------------------------------------------------- Index: /dynamico_lmdz/simple_physics/phyparam/param/comsaison.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/comsaison.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/comsaison.F90 (revision 4176) @@ -0,0 +1,8 @@ +module comsaison + logical,save :: callsais + integer,save :: isaison + real,save :: dist_sol,declin +!$OMP THREADPRIVATE(callsais,isaison,dist_sol,declin) +contains + +end module comsaison Index: /dynamico_lmdz/simple_physics/phyparam/param/constlim.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/constlim.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/constlim.F90 (revision 4176) @@ -0,0 +1,20 @@ +module constlim + REAL,ALLOCATABLE,SAVE :: dgcdrag(:),cdzconst(:) + REAL :: ais1,ais2,cdzmin,ccdzh,cpgam,capsol + REAL :: cdrat,dtradia,lmixmin +!$OMP THREADPRIVATE(dgcdrag,cdzconst) +!$OMP THREADPRIVATE(ais1,ais2,cdzmin,ccdzh,cpgam,capsol) +!$OMP THREADPRIVATE(cdrat,dtradia,lmixmin) +contains + + subroutine InitConstlim + USE mod_phys_lmdz_para + USE dimphy, ONLY : klon,klev + implicit none + + allocate(dgcdrag(klon)) + allocate(cdzconst(klon)) + + end subroutine Initconstlim + +end module constlim Index: /dynamico_lmdz/simple_physics/phyparam/param/convadj.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/convadj.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/convadj.F (revision 4176) @@ -0,0 +1,236 @@ + SUBROUTINE convadj(ngrid,nlay,ptimestep, + S pplay,pplev,ppopsk, + $ pu,pv,ph, + $ pdufi,pdvfi,pdhfi, + $ pduadj,pdvadj,pdhadj) + IMPLICIT NONE + +c======================================================================= +c +c ajustement convectif sec +c on peut ajouter les tendances pdhfi au profil pdh avant l'ajustement +c +c======================================================================= + +c----------------------------------------------------------------------- +c declarations: +c ------------- + +#include "dimensions.h" +#include "comcstfi.h" + +c arguments: +c ---------- + + INTEGER ngrid,nlay + REAL ptimestep + REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay) + REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay) + REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay) + REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay) + +c local: +c ------ + + INTEGER ig,i,l,l1,l2,jj + INTEGER jcnt, jadrs(ngrid) + + REAL*8 sig(nlay+1),sdsig(nlay),dsig(nlay) + REAL*8 zu(ngrid,nlay),zv(ngrid,nlay) + REAL*8 zh(ngrid,nlay) + REAL*8 zu2(ngrid,nlay),zv2(ngrid,nlay) + REAL*8 zh2(ngrid,nlay) + REAL*8 zhm,zsm,zum,zvm,zalpha + + LOGICAL vtest(ngrid),down + +c +c----------------------------------------------------------------------- +c initialisation: +c --------------- +c +c +c----------------------------------------------------------------------- +c detection des profils a modifier: +c --------------------------------- +c si le profil est a modifier +c (i.e. ph(niv_sup) < ph(niv_inf) ) +c alors le tableau "vtest" est mis a .TRUE. ; +c sinon, il reste a sa valeur initiale (.FALSE.) +c cette operation est vectorisable +c On en profite pour copier la valeur initiale de "ph" +c dans le champ de travail "zh" + + + DO 1010 l=1,nlay + DO 1015 ig=1,ngrid + zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep + zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep + zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep +1015 CONTINUE +1010 CONTINUE + + zu2(:,:)=zu(:,:) + zv2(:,:)=zv(:,:) + zh2(:,:)=zh(:,:) + + DO 1020 ig=1,ngrid + vtest(ig)=.FALSE. + 1020 CONTINUE +c + DO 1040 l=2,nlay + DO 1060 ig=1,ngrid +CRAY vtest(ig)=CVMGM(.TRUE. , vtest(ig), +CRAY . zh2(ig,l)-zh2(ig,l-1)) + IF(zh2(ig,l).LT.zh2(ig,l-1)) vtest(ig)=.TRUE. + 1060 CONTINUE + 1040 CONTINUE +c +CRAY CALL WHENNE(ngrid, vtest, 1, 0, jadrs, jcnt) + jcnt=0 + DO 1070 ig=1,ngrid + IF(vtest(ig)) THEN + jcnt=jcnt+1 + jadrs(jcnt)=ig + ENDIF +1070 CONTINUE + + +c----------------------------------------------------------------------- +c Ajustement des "jcnt" profils instables indices par "jadrs": +c ------------------------------------------------------------ +c + DO 1080 jj = 1, jcnt +c + i = jadrs(jj) +c +c Calcul des niveaux sigma sur cette colonne + DO l=1,nlay+1 + sig(l)=pplev(i,l)/pplev(i,1) + ENDDO + DO l=1,nlay + dsig(l)=sig(l)-sig(l+1) + sdsig(l)=ppopsk(i,l)*dsig(l) + ENDDO + l2 = 1 +c +c -- boucle de sondage vers le haut +c +cins$ Loop + 8000 CONTINUE +c + l2 = l2 + 1 +c +cins$ Exit + IF (l2 .GT. nlay) Goto 8001 +c + IF (zh2(i, l2) .LT. zh2(i, l2-1)) THEN +c +c -- l2 est le niveau le plus haut de la colonne instable +c + l1 = l2 - 1 + l = l1 + zsm = sdsig(l2) + zhm = zh2(i, l2) +c +c -- boucle de sondage vers le bas +c +cins$ Loop + 8020 CONTINUE +c + zsm = zsm + sdsig(l) + zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm +c +c -- doit on etendre la colonne vers le bas ? +c +c_EC (M1875) 20/6/87 : AND -> AND THEN +c + down = .FALSE. + IF (l1 .NE. 1) THEN !-- and then + IF (zhm .LT. zh2(i, l1-1)) THEN + down = .TRUE. + END IF + END IF +c + IF (down) THEN +c + l1 = l1 - 1 + l = l1 +c + ELSE +c +c -- peut on etendre la colonne vers le haut ? +c +cins$ Exit + IF (l2 .EQ. nlay) Goto 8021 +c +cins$ Exit + IF (zh2(i, l2+1) .GE. zhm) Goto 8021 +c + l2 = l2 + 1 + l = l2 +c + END IF +c +cins$ End Loop + GO TO 8020 + 8021 CONTINUE +c +c -- nouveau profil : constant (valeur moyenne) +c + zalpha=0. + zum=0. + zvm=0. + DO 1100 l = l1, l2 + zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) + zh2(i, l) = zhm + zum=zum+dsig(l)*zu(i,l) + zvm=zvm+dsig(l)*zv(i,l) + 1100 CONTINUE + zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) + zum=zum/(sig(l1)-sig(l2+1)) + zvm=zvm/(sig(l1)-sig(l2+1)) + IF(zalpha.GT.1.) THEN + PRINT*,'WARNING dans convadj zalpha=',zalpha + if(ig.eq.1) then + print*,'Au pole nord' + elseif (ig.eq.ngrid) then + print*,'Au pole sud' + else + print*,'Point i=', + . ig-((ig-1)/iim)*iim,'j=',(ig-1)/iim+1 + endif + STOP + zalpha=1. + ELSE +c IF(zalpha.LT.0.) STOP'zalpha=0' + IF(zalpha.LT.1.e-5) zalpha=1.e-5 + ENDIF + DO l=l1,l2 + zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) + zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) + ENDDO + + l2 = l2 + 1 +c + END IF +c +cins$ End Loop + GO TO 8000 + 8001 CONTINUE +c + 1080 CONTINUE +c + DO 4000 l=1,nlay + DO 4020 ig=1,ngrid + pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep + pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep + pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep +c pdhadj(ig,l)=0. +c pduadj(ig,l)=0. +c pdvadj(ig,l)=0. + 4020 CONTINUE + 4000 CONTINUE +c + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/ftn00 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/ftn00 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/ftn00 (revision 4176) @@ -0,0 +1,2 @@ + DRS error 22 in routine aslun + Cannot open dictionary file. Index: /dynamico_lmdz/simple_physics/phyparam/param/gather.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/gather.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/gather.F (revision 4176) @@ -0,0 +1,16 @@ + SUBROUTINE monGATHER(n,a,b,index) +c + IMPLICIT NONE +C + INTEGER n,index(n),i + REAL a(n),b(n) +c + DO 100 i=1,n + a(i)=b(index(i)) +100 CONTINUE +c + RETURN + END +c +c voir aussi vec_gather(v,vindices,count,r)...p.11-14 +c Index: /dynamico_lmdz/simple_physics/phyparam/param/inifrict.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/inifrict.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/inifrict.F (revision 4176) @@ -0,0 +1,66 @@ + SUBROUTINE inifrict(timestep) + USE comgeomfi, ONLY : ngridmax,nlayermx + USE constlim + IMPLICIT NONE +c======================================================================= +c +c Calcul des coefficients de ls diffusion verticale +c +c======================================================================= +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "comcstfi.h" +! include "comvert.h" + +c local: +c ------ + + REAL dgrad,cpr,rl2,rrr,timestep + INTEGER l,ij + +c----------------------------------------------------------------------- + + call initconstlim + + dtradia=timestep + PRINT*,'DTPHYS',dtradia + lmixmin=100. + ais1 = 1. + ais2 = ais1 - 1. + print*,'ais1',ais1,'ais2',ais2 + cdzmin = 1.e-6 + OPEN(99,file='cdzmin',status='old',err=9999) + READ(99,*) cdzmin +9999 CLOSE(99) + PRINT*,'cdzmin=',cdzmin + + cpr = cpp/ r + ccdzh = 2.5*g +c!!! cpgam = 5.e-3*cpp + cpgam=0. + +c----------------------------------------------------------------------- +c coefficient de diffusion dans l'atmosphere: +c ------------------------------------------- + + rl2=lmixmin**2 + cdzconst(1)= 0. + DO 15 l=1,nlayermx - 1 + cdzconst(l+1)= dtradia*g*g*cpr*rl2 + print*,'cdzconst(',l+1,') = ',cdzconst(l+1) + 15 CONTINUE + +c----------------------------------------------------------------------- +c couche limite de surface: +c ------------------------- + + cdrat = 2.e-3 + dgrad = dtradia*g*cpp/r + DO 16 ij = 1, ngridmax + dgcdrag( ij ) = cdrat * dgrad + 16 CONTINUE + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/iniorbit.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iniorbit.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iniorbit.F (revision 4176) @@ -0,0 +1,100 @@ + SUBROUTINE iniorbit + $ (paphelie,pperiheli,pyear_day,pperi_day,pobliq) + IMPLICIT NONE + +c======================================================================= +c +c Auteur: +c ------- +c Frederic Hourdin 22 Fevrier 1991 +c +c Objet: +c ------ +c Initialisation du sous programme orbite qui calcule +c a une date donnee de l'annee de duree year_day commencant +c a l'equinoxe de printemps et dont le perihelie se situe +c a la date peri_day, la distance au soleil et la declinaison. +c +c Interface: +c ---------- +c - Doit etre appele avant d'utiliser orbite. +c - initialise le common comorbit +c +c Arguments: +c ---------- +c +c Input: +c ------ +c aphelie \ aphelie et perihelie de l'orbite +c periheli / en millions de kilometres. +c +c======================================================================= + +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "comorbit.h" + +c Arguments: +c ---------- + + REAL paphelie,pperiheli,pyear_day,pperi_day,pobliq + +c Local: +c ------ + + REAL zxref,zanom,zz,zx0,zdx + INTEGER iter + +c----------------------------------------------------------------------- + + aphelie =paphelie + periheli=pperiheli + year_day=pyear_day + obliquit=pobliq + peri_day=pperi_day + + pi=2.*asin(1.) + PRINT*,'Perihelie en Mkm ',periheli + PRINT*,'Aphelise en Mkm ',aphelie + PRINT*,'obliquite en degres :',obliquit + unitastr=149.597927 + e_elips=(aphelie-periheli)/(periheli+aphelie) + p_elips=0.5*(periheli+aphelie)*(1-e_elips*e_elips)/unitastr + + print*,'e_elips',e_elips + print*,'p_elips',p_elips + +c----------------------------------------------------------------------- +c calcul de l'angle polaire et de la distance au soleil : +c ------------------------------------------------------- + +c calcul de l'zanomalie moyenne + + zz=(year_day-pperi_day)/year_day + zanom=2.*pi*(zz-nint(zz)) + zxref=abs(zanom) + PRINT*,'zanom ',zanom + +c resolution de l'equation horaire zx0 - e * sin (zx0) = zxref +c methode de Newton + + zx0=zxref+e_elips*sin(zxref) + DO 110 iter=1,100 + zdx=-(zx0-e_elips*sin(zx0)-zxref)/(1.-e_elips*cos(zx0)) + if(abs(zdx).le.(1.e-12)) goto 120 + zx0=zx0+zdx +110 continue +120 continue + zx0=zx0+zdx + if(zanom.lt.0.) zx0=-zx0 + PRINT*,'zx0 ',zx0 + +c zteta est la longitude solaire + + timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) + PRINT*,'longitude solaire du perihelie timeperi = ',timeperi + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F (revision 4176) @@ -0,0 +1,211 @@ + SUBROUTINE iniphyparam(ngrid,nlayer, + $ punjours, + $ pdayref,ptimestep, + $ prad,pg,pr,pcpp) + use IOIPSL + use getparam + use dimphy + USE mod_grid_phy_lmdz + USE mod_phys_lmdz_para + use comgeomfi + use comsaison + USE geometry_mod, ONLY : longitude,latitude,cell_area + + IMPLICIT NONE + +c +c======================================================================= +c +c subject: +c -------- +c +c Initialisation for the physical parametrisations of the LMD +c martian atmospheric general circulation modele. +c +c author: Frederic Hourdin 15 / 10 /93 +c ------- +c +c arguments: +c ---------- +c +c input: +c ------ +c +c ngrid Size of the horizontal grid. +c All internal loops are performed on that grid. +c nlayer Number of vertical layers. +c pdayref Day of reference for the simulation +c firstcall True at the first call +c lastcall True at the last call +c pday Number of days counted from the North. Spring +c equinoxe. +c +c======================================================================= +c +c----------------------------------------------------------------------- +c declarations: +c ------------- + + +#include "planete.h" +#include "comcstfi.h" +#include "callkeys.h" +#include "surface.h" +#include "iniprint.h" + + + REAL prad,pg,pr,pcpp,punjours + + INTEGER ngrid,nlayer + REAL pdayref + + REAL ptimestep + INTEGER ig,ierr,offset + + EXTERNAL inifrict,iniorbit,orbite + + print*,'INIPHYPARAM' + CALL InitComgeomfi + + IF (klon.NE.klon_omp) THEN + PRINT*,'STOP in iniphyparam' + PRINT*,'Probleme de dimenesions :' + PRINT*,'klon = ',klon + PRINT*,'klon_omp = ',klon_omp + STOP + ENDIF + + IF (nlayer.NE.nlayermx) THEN + PRINT*,'STOP in iniphyparam' + PRINT*,'Probleme de dimenesions :' + PRINT*,'nlayer = ',nlayer + PRINT*,'nlayermx = ',nlayermx + STOP + ENDIF + + IF (ngrid.NE.klon_glo) THEN + PRINT*,'STOP in iniphyparam' + PRINT*,'Probleme de dimenesions :' + PRINT*,'ngrid = ',ngrid + PRINT*,'ngridmax = ',klon_glo +! STOP + ENDIF + + print*,'Avant les getpar ' + CALL getpar('unjours',86400. ,unjours,'unjours') + CALL getpar('rad',6400000. ,rad,'rad') + CALL getpar('g',9.8 ,g,'g') + CALL getpar('cpp',1004. ,cpp,'cpp') + CALL getpar('mugaz',28. ,mugaz,'mugaz') + CALL getpar('year_day',360. ,year_day,'year_day') + CALL getpar('periheli',150. ,periheli,'periheli') + CALL getpar('aphelie',150. ,aphelie,'aphelie') + CALL getpar('peri_day',0. ,peri_day,'peri_day') + CALL getpar('obliquit',23. ,obliquit,'obliquit') + CALL getpar('Cd_mer',.01 ,Cd_mer,'Cd_mer') + CALL getpar('Cd_ter',.01 ,Cd_ter,'Cd_ter') + CALL getpar('I_mer',30000. ,I_mer,'I_mer') + CALL getpar('I_ter',30000. ,I_ter,'I_ter') + CALL getpar('alb_ter',.112 ,alb_ter,'alb_ter') + CALL getpar('alb_mer',.112 ,alb_mer,'alb_mer') + CALL getpar('emi_mer',1. ,emi_mer,'emi_mer') + CALL getpar('emi_mer',1. ,emi_mer,'emi_mer') + CALL getpar('emin_turb',1.e-16 ,emin_turb,'emin_turb') + CALL getpar('lmixmin',100. ,lmixmin,'lmixmin') + CALL getpar('coefvis',.99 ,coefvis,'coefvis') + CALL getpar('coefir',.08 ,coefir,'coefir') + + + CALL getpar('callrad',.true.,callrad,'appel rayonnemen') + CALL getpar('calldifv',.true.,calldifv,'appel difv') + CALL getpar('calladj',.true.,calladj,'appel adj') + CALL getpar('callcond',.true.,callcond,'appel cond') + CALL getpar('callsoil',.true.,callsoil,'appel soil') + CALL getpar('season',.true.,season,'appel soil') + CALL getpar('diurnal',.false.,diurnal,'appel soil') + CALL getpar('lverbose',.true.,lverbose,'appel soil') + CALL getpar('period_sort',1.,period_sort,'period sorties en jour') + + write(lunout,*) 'unjours=',unjours + write(lunout,*) 'rad=',rad + write(lunout,*) 'g=',g + write(lunout,*) 'cpp=',cpp + write(lunout,*) 'mugaz=',mugaz + write(lunout,*) 'year_day=',year_day + write(lunout,*) 'periheli=',periheli + write(lunout,*) 'aphelie=',aphelie + write(lunout,*) 'peri_day=',peri_day + write(lunout,*) 'obliquit=',obliquit + write(lunout,*) 'Cd_mer=',Cd_mer + write(lunout,*) 'Cd_ter=',Cd_ter + write(lunout,*) 'I_mer=',I_mer + write(lunout,*) 'I_ter=',I_ter + write(lunout,*) 'alb_ter=',alb_ter + write(lunout,*) 'alb_mer=',alb_mer + write(lunout,*) 'emi_mer=',emi_mer + write(lunout,*) 'emi_mer=',emi_mer + write(lunout,*) 'emin_turb=',emin_turb + write(lunout,*) 'lmixmin=',lmixmin + write(lunout,*) 'coefvis=',coefvis + write(lunout,*) 'coefir=',coefir + write(lunout,*) 'callrad=',callrad + write(lunout,*) 'calldifv=',calldifv + write(lunout,*) 'calladj=',calladj + write(lunout,*) 'callcond=',callcond + write(lunout,*) 'callsoil=',callsoil + write(lunout,*) 'season=',season + write(lunout,*) 'diurnal=',diurnal + write(lunout,*) 'lverbose=',lverbose + write(lunout,*) 'period_sort=',period_sort + + print*,'Activation de la physique:' + print*,' Rayonnement ',callrad + print*,' Diffusion verticale turbulente ', calldifv + print*,' Ajustement convectif ',calladj + print*,' Sol ',callsoil + print*,' Cycle diurne ',diurnal + +c choice of the frequency of the computation of radiations + IF(diurnal) THEN + iradia=NINT(punjours/(20.*ptimestep)) + ELSE + iradia=NINT(punjours/(4.*ptimestep)) + ENDIF + iradia=1 + PRINT*,'unjours',punjours + PRINT*,'The radiative transfer is computed each ', + s iradia,' physical time-step or each ', + s iradia*ptimestep,' seconds' +c----------------------------------------------------------------------- + + offset=klon_mpi_begin-1 + + print*,'latitude0 ohe',latitude(1:3),latitude(klon_omp) +! long(1:klon_omp)=plon(offset+klon_omp_begin:offset+klon_omp_end) +! lati(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end) +! area(1:klon_omp)=parea(offset+klon_omp_begin:offset+klon_omp_end) + long(1:klon_omp)=longitude(1:klon_omp) + lati(1:klon_omp)=latitude(1:klon_omp) + area(1:klon_omp)=cell_area(1:klon_omp) + totarea=sum(cell_area,ngrid) + print*,'OK17 AAA' + + sinlat(:)=sin(lati(:)) + coslat(:)=cos(lati(:)) + sinlon(:)=sin(long(:)) + coslon(:)=cos(long(:)) + + pi=2.*asin(1.) + prad=rad + pg=g + r=8.134/(mugaz*0.001) + print*,'R=',r + pr=r + pcpp=cpp + rcp=r/cpp + dtphys=ptimestep + punjours=unjours + + RETURN +9999 STOP'Cette version demande les fichier rnatur.dat et surf.def' + END Index: /dynamico_lmdz/simple_physics/phyparam/param/iniphysiq_param.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iniphysiq_param.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iniphysiq_param.F (revision 4176) @@ -0,0 +1,129 @@ +! +! $Id: iniphysiq.F 1403 2010-07-01 09:02:53Z fairhead $ +! +c +c + SUBROUTINE iniphysiq_param(ngrid,nlayer, + $ punjours, + $ pdayref,ptimestep, + $ plat,plon,parea,pcu,pcv, + $ prad,pg,pr,pcpp,iflag_phys) + USE dimphy + USE mod_grid_phy_lmdz + USE mod_phys_lmdz_para + USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg +! USE inigeomphy_mod, ONLY : lonfi,latfi + + + IMPLICIT NONE +c +c======================================================================= +c +c subject: +c -------- +c +c Initialisation for the physical parametrisations of the LMD +c martian atmospheric general circulation modele. +c +c author: Frederic Hourdin 15 / 10 /93 +c ------- +c +c arguments: +c ---------- +c +c input: +c ------ +c +c ngrid Size of the horizontal grid. +c All internal loops are performed on that grid. +c nlayer Number of vertical layers. +c pdayref Day of reference for the simulation +c firstcall True at the first call +c lastcall True at the last call +c pday Number of days counted from the North. Spring +c equinoxe. +c +c======================================================================= +c +c----------------------------------------------------------------------- +c declarations: +c ------------- + +#include "comcstphy.h" + REAL prad,pg,pr,pcpp,punjours + + INTEGER ngrid,nlayer,iflag_phys + REAL plat(ngrid),plon(ngrid),parea(klon_glo) + REAL pcu(klon_glo),pcv(klon_glo) + INTEGER pdayref + INTEGER :: ibegin,iend,offset,indmin,indmax + REAL pi + + REAL ptimestep + CHARACTER (LEN=20) :: modname='iniphysiq' + CHARACTER (LEN=80) :: abort_message + + pi=2.*asin(1.) + print*,'INInnn ENTREE DANS INIPHYSIQ' + IF (nlayer.NE.klev) THEN + PRINT*,'STOP in inifis' + PRINT*,'Probleme de dimensions :' + PRINT*,'nlayer = ',nlayer + PRINT*,'klev = ',klev + abort_message = '' + CALL abort_gcm (modname,abort_message,1) + ENDIF + + IF (ngrid.NE.klon_omp) THEN + PRINT*,'STOP in inifis' + PRINT*,'Probleme de dimensions :' + PRINT*,'ngrid = ',ngrid + PRINT*,'klon = ',klon_omp + abort_message = '' + CALL abort_gcm (modname,abort_message,1) + ENDIF +c$OMP PARALLEL PRIVATE(ibegin,iend) +c$OMP+ SHARED(parea,pcu,pcv,plon,plat) + + print*,'Dans iniphysiq ' + offset=klon_mpi_begin-1 +! cell_area(1:klon_omp)=parea(offset+klon_omp_begin: +! & offset+klon_omp_end) + !cuphy(1:klon_omp)=pcu(offset+klon_omp_begin:offset+klon_omp_end) + !cvphy(1:klon_omp)=pcv(offset+klon_omp_begin:offset+klon_omp_end) + indmin=offset+klon_omp_begin + indmax=offset+klon_omp_end +! longitude_deg(1:klon_omp)=180./pi*lonfi(indmin:indmax) +! latitude_deg(1:klon_omp)=180./pi*latfi(indmin:indmax) + print*,'latitude0 deg ',latitude_deg(1),latitude_deg(klon_omp) + +! call suphel +! prad,pg,pr,pcpp + rradius=prad + rg=pg + rr=pr + rcpp=pcpp + +! return + + CALL iniphyparam(ngrid,nlayer, + $ punjours, + $ pdayref,ptimestep, + $ prad,pg,pr,pcpp) + + + print*,'OK iniphyparam ',size(plat) + +c$OMP END PARALLEL + + print*,'ATTENTION !!! TRAVAILLER SUR INIPHYSIQ' + print*,'CONTROLE DES LATITUDES, LONGITUDES, PARAMETRES ...' + + RETURN +9999 CONTINUE + abort_message ='Cette version demande les fichier rnatur.dat + & et surf.def' + CALL abort_gcm (modname,abort_message,1) + + + END Index: /dynamico_lmdz/simple_physics/phyparam/param/interpol.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/interpol.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/interpol.F (revision 4176) @@ -0,0 +1,95 @@ + SUBROUTINE interpol(xinp,yinp,ninp,xout,yout,nout) + IMPLICIT NONE + +c======================================================================= +c +c objet: +c interpolation lineaire +c +c ATTENTION!: +c les tableaux x doivent etre croissants +c +c Arguments: +c ---------- +c +c entree: +c ------- +c ninp INT nombre de points des donnes +c xinp(ninp) REAL abbcisses des donnees (doivent etre croissantes) +c yinp(ninp) REAL valeurs des donnees +c +c sorties: +c -------- +c nout INT nombre de points interpolles +c xout(nout) REAL abbcisses (doivent etre croissantes) +c yout(nout) REAL valeurs interpollees +c +c======================================================================= + +c----------------------------------------------------------------------- +c 0. Declarations : +c ----------------- + + INTEGER ninp,nout + INTEGER iinp,iout,iout0,iout1 + REAL xinp(ninp),yinp(ninp) + REAL xout(nout),yout(nout) + +c----------------------------------------------------------------------- +c traitement des valeurs telles que xoutxinp(ninp): +c ----------------------------------------------- + + DO 500 iout=iout1,nout + yout(iout)=yinp(ninp) +c si on veut extrapoler aux bords: +c 1 +(xout(iout)-xinp(ninp))* +c 2 (yinp(ninp-1)-yinp(ninp))/(xinp(ninp-1)-xinp(ninp)) +500 CONTINUE + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/iophys.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iophys.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iophys.F90 (revision 4176) @@ -0,0 +1,110 @@ + subroutine iophys_ecrit(nom,lllm,titre,unite,px) + USE dimphy + USE mod_phys_lmdz_para + USE mod_grid_phy_lmdz + IMPLICIT NONE + + + +! Ecriture de variables diagnostiques au choix dans la physique +! dans un fichier NetCDF nomme 'diagfi'. Ces variables peuvent etre +! 3d (ex : temperature), 2d (ex : temperature de surface), ou +! 0d (pour un scalaire qui ne depend que du temps : ex : la longitude +! solaire) +! (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne) +! La periode d'ecriture est donnee par +! "ecritphy " regle dans le fichier de controle de run : run.def +! +! writediagfi peut etre appele de n'importe quelle subroutine +! de la physique, plusieurs fois. L'initialisation et la creation du +! fichier se fait au tout premier appel. +! +! WARNING : les variables dynamique (u,v,t,q,ps) +! sauvees par writediagfi avec une +! date donnee sont legerement differentes que dans le fichier histoire car +! on ne leur a pas encore ajoute de la dissipation et de la physique !!! +! IL est RECOMMANDE d'ajouter les tendance physique a ces variables +! avant l'ecriture dans diagfi (cf. physiq.F) +! +! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4 +! +! parametres (input) : +! ---------- +! unit : unite logique du fichier de sortie (toujours la meme) +! nom : nom de la variable a sortir (chaine de caracteres) +! titre: titre de la variable (chaine de caracteres) +! unite : unite de la variable (chaine de caracteres) +! px : variable a sortir (real 0, 1, 2, ou 3d) +! +!================================================================= + +#include "dimensions.h" +#include "paramet.h" +#include "netcdf.inc" +#include "iotd.h" + + +! Arguments on input: + integer lllm + character (len=*) :: nom,titre,unite + integer imjmax + parameter (imjmax=100000) + real px(klon_omp,lllm) + real xglo(klon_glo,lllm) + real zx(iim,jjp1,lllm) + + + CALL Gather(px,xglo) +!$OMP MASTER + IF (is_mpi_root) THEN + CALL Grid1Dto2D_glo(xglo,zx) + call iotd_ecrit(nom,lllm,titre,unite,zx) + ENDIF +!$OMP END MASTER + + return + end + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + SUBROUTINE iophys_ini(fichnom,presnivs) + USE mod_phys_lmdz_para + + IMPLICIT NONE + +!======================================================================= +! +! Auteur: L. Fairhead , P. Le Van, Y. Wanherdrick, F. Forget +! ------- +! +! Objet: +! ------ +! +! 'Initialize' the diagfi.nc file: write down dimensions as well +! as time-independent fields (e.g: geopotential, mesh area, ...) +! +!======================================================================= +!----------------------------------------------------------------------- +! Declarations: +! ------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comgeom.h" +! #include "comvert.h" +REAL presnivs(llm) + +character*10 fichnom +real pi + +! Arguments: +! ---------- + +!$OMP MASTER + IF (is_mpi_root) THEN +pi=2.*asin(1.) +call iotd_ini('phys.nc ', & +iim,jjp1,llm,rlonv(1:iim)*180./pi,rlatu*180./pi,presnivs) + ENDIF +!$OMP END MASTER + + END Index: /dynamico_lmdz/simple_physics/phyparam/param/iotd.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iotd.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iotd.h (revision 4176) @@ -0,0 +1,16 @@ +!======================================================================= +! +! Auteur: F. Hourdin +! ------- +! +! Objet: +! ------ +! Light interface for netcdf outputs. can be used outside LMDZ +! +!======================================================================= + + integer imax,jmax,lmax,nid + INTEGER dim_coord(4) + real iotd_ts + + common/ecritd_c/imax,jmax,lmax,nid,dim_coord,iotd_ts Index: /dynamico_lmdz/simple_physics/phyparam/param/iotd_ecrit.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iotd_ecrit.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iotd_ecrit.F90 (revision 4176) @@ -0,0 +1,180 @@ + subroutine iotd_ecrit(nom,llm,titre,unite,px) + + +!======================================================================= +! +! Auteur: F. Hourdin +! ------- +! +! Objet: +! ------ +! Light interface for netcdf outputs. can be used outside LMDZ +! +!======================================================================= +!----------------------------------------------------------------------- +! ---------- +! nom : nom de la variable a sortir (chaine de caracteres) +! llm : nombre de couches +! titre: titre de la variable (chaine de caracteres) +! unite : unite de la variable (chaine de caracteres) +! px : variable a sortir +! +!================================================================= + + implicit none + +! Commons + +#include "netcdf.inc" +#include "iotd.h" + + +! Arguments on input: + integer llm + character (len=*) :: nom,titre,unite + integer imjmax + parameter (imjmax=100000) + real px(imjmax*llm) + +! Local variables: + + real*4 date + real*4 zx(imjmax*llm) + + + integer ierr,ndim,dim_cc(4) + integer iq + integer i,j,l + + integer zitau + character firstnom*20 + SAVE firstnom + SAVE zitau + SAVE date + data firstnom /'1234567890'/ + data zitau /0/ + +! Ajouts + integer, save :: ntime=0 + integer :: idim,varid + character (len =50):: fichnom + integer, dimension(4) :: id + integer, dimension(4) :: edges,corner + + +!*************************************************************** +! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file +! ------------------------------------------------------------------------ +! (Au tout premier appel de la subroutine durant le run.) + + +!-------------------------------------------------------- +! Write the variables to output file if it's time to do so +!-------------------------------------------------------- + + +! Compute/write/extend 'Time' coordinate (date given in days) +! (done every "first call" (at given time level) to writediagfi) +! Note: date is incremented as 1 step ahead of physics time +!-------------------------------------------------------- + + zx(1:imax*jmax*llm)=px(1:imax*jmax*llm) + if (firstnom =='1234567890') then + firstnom=nom + endif + +! print*,'nom ',nom,firstnom + +!! Quand on tombe sur la premiere variable on ajoute un pas de temps + if (nom.eq.firstnom) then + ! We have identified a "first call" (at given date) + + ntime=ntime+1 ! increment # of stored time steps + +!! print*,'ntime ',ntime + date=ntime +! date= float (zitau +1)/float (day_step) + + ! compute corresponding date (in days and fractions thereof) + ! Get NetCDF ID of 'Time' variable + + ierr=NF_SYNC(nid) + + ierr= NF_INQ_VARID(nid,"Time",varid) + ! Write (append) the new date to the 'Time' array + + + ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date) + +! print*,'date ',date,ierr,nid +! print*,'IOTD Date ,varid,nid,ntime,date',varid,nid,ntime,date + + if (ierr.ne.NF_NOERR) then + write(*,*) "***** PUT_VAR matter in writediagfi_nc" + write(*,*) "***** with time" + write(*,*) 'ierr=', ierr + endif + +! write(6,*)'WRITEDIAGFI: date= ', date + end if ! of if (nom.eq.firstnom) + + +!Case of a 3D variable +!--------------------- + if (llm==lmax) then + ndim=4 + corner(1)=1 + corner(2)=1 + corner(3)=1 + corner(4)=ntime + edges(1)=imax + edges(2)=jmax + edges(3)=llm + edges(4)=1 + dim_cc=dim_coord + + +!Case of a 2D variable +!--------------------- + + else if (llm==1) then + ndim=3 + corner(1)=1 + corner(2)=1 + corner(3)=ntime + corner(4)=1 + edges(1)=imax + edges(2)=jmax + edges(3)=1 + edges(4)=1 + dim_cc(1:2)=dim_coord(1:2) + dim_cc(3)=dim_coord(4) + + endif ! of if llm=1 ou llm + +! AU premier pas de temps, on crée les variables +!----------------------------------------------- + + if (ntime==1) then + ierr = NF_REDEF (nid) + ierr = NF_DEF_VAR(nid,nom,NF_FLOAT,ndim,dim_cc,varid) + print*,'DEF ',nom,nid,varid + ierr = NF_ENDDEF(nid) + else + ierr= NF_INQ_VARID(nid,nom,varid) + print*,'INQ ',nom,nid,varid +! Commandes pour recuperer automatiquement les coordonnees +! ierr= NF_INQ_DIMID(nid,"longitude",id(1)) + endif + + + ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,zx) + + if (ierr.ne.NF_NOERR) then + write(*,*) "***** PUT_VAR problem in writediagfi" + write(*,*) "***** with ",nom + write(*,*) 'ierr=', ierr + endif + + + end Index: /dynamico_lmdz/simple_physics/phyparam/param/iotd_fin.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iotd_fin.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iotd_fin.F90 (revision 4176) @@ -0,0 +1,25 @@ + SUBROUTINE iotd_fin + IMPLICIT NONE + +!======================================================================= +! +! Auteur: F. Hourdin +! ------- +! +! Objet: +! ------ +! Light interface for netcdf outputs. can be used outside LMDZ +! +!======================================================================= + + +#include "netcdf.inc" +#include "iotd.h" + integer ierr + +! Arguments: +! ---------- + + ierr=NF_close(nid) + + END Index: /dynamico_lmdz/simple_physics/phyparam/param/iotd_ini.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/iotd_ini.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/iotd_ini.F90 (revision 4176) @@ -0,0 +1,116 @@ + SUBROUTINE iotd_ini(fichnom,iim,jjm,llm,prlonv,prlatu,pcoordv) + IMPLICIT NONE + +!======================================================================= +! +! Auteur: F. Hourdin +! ------- +! +! Objet: +! ------ +! Light interface for netcdf outputs. can be used outside LMDZ +! +!======================================================================= +!----------------------------------------------------------------------- +! Declarations: +! ------------- + +#include "netcdf.inc" +#include "iotd.h" + +! Arguments: +! ---------- + + integer iim,jjm,llm + real prlonv(iim),prlatu(jjm),pcoordv(llm),timestep + INTEGER id_FOCE + + integer corner(4),edges(4),ndim + real px(1000) + character (len=10) :: nom + +! Local: +! ------ + INTEGER ierr + + integer :: nvarid + integer, dimension(2) :: id + integer :: varid + + character*10 fichnom + real*4 rlonv(iim),rlatu(jjm),coordv(llm) + + real pi + + print*,'INIIO prlonv ',prlonv + imax=iim + jmax=jjm + lmax=llm + + rlonv=prlonv + rlatu=prlatu + coordv=pcoordv + +!----------------------------------------------------------------------- + + + pi=2.*asin(1.) + +! Define dimensions + + ! Create the NetCDF file + ierr=NF_CREATE(fichnom, NF_CLOBBER, nid) + ! Define the 'Time' dimension + ierr=nf_def_dim(nid,"Time",NF_UNLIMITED,dim_coord(4)) + ! Define the 'Time' variable + ierr=NF_DEF_VAR(nid, "Time", NF_FLOAT, 1, dim_coord(4),varid) +! ! Add a long_name attribute +! ierr=NF_PUT_ATT_TEXT(nid, varid, "long_name",4,"Time") +! ! Add a units attribute + ierr=NF_PUT_ATT_TEXT(nid, varid,'units',29,"days since 0000-00-0 00:00:00") + ! Switch out of NetCDF Define mode + + ierr=NF_DEF_DIM(nid, "longitude", iim, dim_coord(1)) + ierr=NF_DEF_DIM(nid, "latitude", jjm, dim_coord(2)) + ierr=NF_DEF_DIM(nid, "altitude", llm, dim_coord(3)) + + + ierr=NF_ENDDEF(nid) +! +! Contol parameters for this run +! -------------------------- + + ierr=NF_REDEF(nid) + ierr=NF_DEF_VAR(nid,"longitude", NF_FLOAT, 1, dim_coord(1),nvarid) +! ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name", 14, +! . "East longitude") +! ierr=NF_PUT_ATT_TEXT(nid,nvarid,'units',12,"degrees_east") + ierr=NF_ENDDEF(nid) + ierr=NF_PUT_VAR_REAL(nid,nvarid,rlonv) + print*,ierr + +! -------------------------- + ierr=NF_REDEF(nid) + ierr=NF_DEF_VAR(nid, "latitude", NF_FLOAT, 1, dim_coord(2),nvarid) +! ierr=NF_PUT_ATT_TEXT(nid,nvarid,'units',13,"degrees_north") +! ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name", 14,"North latitude") + ierr=NF_ENDDEF(nid) + ierr=NF_PUT_VAR_REAL(nid,nvarid,rlatu) +! +! -------------------------- + ierr=NF_REDEF(nid) + ierr=NF_DEF_VAR(nid, "altitude", NF_FLOAT, 1,dim_coord(3),nvarid) + ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",10,"pseudo-alt") +! ierr=NF_PUT_ATT_TEXT(nid,nvarid,'units',2,"km") + if ( pcoordv(2)>pcoordv(1) ) then + ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",10,"pseudo-alt") + ierr=NF_PUT_ATT_TEXT(nid,nvarid,'positive',2,"up") + else + ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",8,"pressure") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',4,"down") + endif + ierr=NF_ENDDEF(nid) + + ierr=NF_PUT_VAR_REAL(nid,nvarid,coordv) +! + END Index: /dynamico_lmdz/simple_physics/phyparam/param/lw.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/lw.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/lw.F (revision 4176) @@ -0,0 +1,220 @@ + SUBROUTINE lw(ngrid,nlayer,coefir,emissiv, + $ pp,ps_rad,ptsurf,pt, + $ pfluxir,pdtlw, + $ lwrite) + + IMPLICIT NONE +c======================================================================= +c +c calcul de l'evolution de la temperature sous l'effet du rayonnement +c infra-rouge. +c Pour simplifier, les transmissions sont precalculees et ne +c dependent que de l'altitude. +c +c arguments: +c ---------- +c +c entree: +c ------- +c ngrid nombres de points de la grille horizontale +c nlayer nombre de couches +c ptsurf(ngrid) temperature de la surface +c pt(ngrid,nlayer) temperature des couches +c pp(ngrid,nlayer+1) pression entre les couches +c lwrite variable logique pour sorties +c +c sortie: +c ------- +c pdtlw(ngrid,nlayer) taux de refroidissement +c pfluxir(ngrid) flux infrarouge sur le sol +c +c======================================================================= + +c declarations: +c ------------- + +#include "comcstfi.h" + +c arguments: +c ---------- + + INTEGER ngrid,nlayer + REAL coefir,emissiv(ngrid),ps_rad + REAL ptsurf(ngrid),pt(ngrid,nlayer),pp(ngrid,nlayer+1) + REAL pdtlw(ngrid,nlayer),pfluxir(ngrid) + LOGICAL lwrite + +c variables locales: +c ------------------ + + INTEGER nlevel,ilev,ig,i,il + REAL zplanck(ngrid,nlayer+1),zcoef + REAL zfluxup(ngrid,nlayer+1),zfluxdn(ngrid,nlayer+1) + REAL zflux(ngrid,nlayer+1) + REAL zlwtr1(ngrid),zlwtr2(ngrid) + REAL zup(ngrid,nlayer+1),zdup(ngrid) + REAL stephan + + LOGICAL lstrong + SAVE lstrong,stephan + DATA lstrong/.true./ +!$OMP THREADPRIVATE(lstrong,stephan) + +c----------------------------------------------------------------------- +c initialisations: +c ---------------- + + nlevel=nlayer+1 + stephan=5.67e-08 + +c----------------------------------------------------------------------- +c 2. calcul des quantites d'absorbants: +c ------------------------------------- + +c absorption forte + IF(lstrong) THEN + DO ilev=1,nlevel + DO ig=1,ngrid + zup(ig,ilev)=pp(ig,ilev)*pp(ig,ilev)/(2.*g) + ENDDO + ENDDO + IF(lwrite) THEN + DO ilev=1,nlayer + PRINT*,' up(',ilev,') = ',zup(ngrid/2+1,ilev) + ENDDO + ENDIF + zcoef=-log(coefir)/sqrt(ps_rad*ps_rad/(2.*g)) + +c absorption faible + ELSE + DO ilev=1,nlevel + DO ig=1,ngrid + zup(ig,ilev)=pp(ig,ilev) + ENDDO + ENDDO + zcoef=-log(coefir)/ps_rad + ENDIF + + +c----------------------------------------------------------------------- +c 2. calcul de la fonction de corps noir: +c --------------------------------------- + + DO ilev=1,nlayer + DO ig=1,ngrid + zplanck(ig,ilev)=pt(ig,ilev)*pt(ig,ilev) + zplanck(ig,ilev)=stephan* + $ zplanck(ig,ilev)*zplanck(ig,ilev) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 4. flux descendants: +c -------------------- + + DO ilev=1,nlayer + DO ig=1,ngrid + zfluxdn(ig,ilev)=0. + ENDDO + DO ig=1,ngrid + zdup(ig)=zup(ig,ilev)-zup(ig,nlevel) + ENDDO + CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1) + + DO il=nlayer,ilev,-1 + zlwtr2(:)=zlwtr1(:) + DO ig=1,ngrid + zdup(ig)=zup(ig,ilev)-zup(ig,il) + ENDDO + CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1) + DO ig=1,ngrid + zfluxdn(ig,ilev)=zfluxdn(ig,ilev)+ + $ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig)) + ENDDO + ENDDO + ENDDO + + DO ig=1,ngrid + zfluxdn(ig,nlevel)=0. + pfluxir(ig)=emissiv(ig)*zfluxdn(ig,1) + ENDDO + + DO ig=1,ngrid + zfluxup(ig,1)=ptsurf(ig)*ptsurf(ig) + zfluxup(ig,1)=emissiv(ig)*stephan*zfluxup(ig,1)*zfluxup(ig,1) + $ +(1.-emissiv(ig))*zfluxdn(ig,1) + ENDDO + +c----------------------------------------------------------------------- +c 3. flux montants: +c ------------------ + + DO ilev=1,nlayer + DO ig=1,ngrid + zdup(ig)=zup(ig,1)-zup(ig,ilev+1) + ENDDO + CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1) + DO ig=1,ngrid + zfluxup(ig,ilev+1)=zfluxup(ig,1)*zlwtr1(ig) + ENDDO + DO il=1,ilev + zlwtr2(:)=zlwtr1(:) + DO ig=1,ngrid + zdup(ig)=zup(ig,il+1)-zup(ig,ilev+1) + ENDDO + CALL lwtr(ngrid,zcoef,lstrong,zdup,zlwtr1) + DO ig=1,ngrid + zfluxup(ig,ilev+1)=zfluxup(ig,ilev+1)+ + $ zplanck(ig,il)*(zlwtr1(ig)-zlwtr2(ig)) + ENDDO + ENDDO + + ENDDO + +c----------------------------------------------------------------------- +c 5. calcul des flux nets: +c ------------------------ + + DO ilev=1,nlevel + DO ig=1,ngrid + zflux(ig,ilev)=zfluxup(ig,ilev)-zfluxdn(ig,ilev) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 6. Calcul des taux de refroidissement: +c -------------------------------------- + + DO ilev=1,nlayer + DO ig=1,ngrid + pdtlw(ig,ilev)=(zflux(ig,ilev+1)-zflux(ig,ilev))* + $ g/(cpp*(pp(ig,ilev+1)-pp(ig,ilev))) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 10. sorties eventuelles: +c ------------------------ + + IF (lwrite) THEN + + PRINT* + PRINT*,'Diagnostique rayonnement thermique' + PRINT* + PRINT*,'temperature ', + $ 'flux montant flux desc. taux de refroid.' + i=ngrid/2+1 + WRITE(6,9000) ptsurf(i) + DO ilev=1,nlayer + WRITE(6,'(i4,4e18.4)') ilev,pt(i,ilev), + $ zfluxup(i,ilev),zfluxdn(i,ilev),pdtlw(i,ilev) + ENDDO + WRITE(6,9000) zfluxup(i,nlevel),zfluxdn(i,nlevel) + + ENDIF + +c----------------------------------------------------------------------- + + RETURN +9000 FORMAT(4e18.4) + END Index: /dynamico_lmdz/simple_physics/phyparam/param/lwtr.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/lwtr.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/lwtr.F (revision 4176) @@ -0,0 +1,18 @@ + SUBROUTINE lwtr(ngrid,coef,lstrong,dup,transm) + IMPLICIT NONE + INTEGER ngrid + REAL coef + LOGICAL lstrong + REAL dup(ngrid),transm(ngrid) + INTEGER ig + IF(lstrong) THEN + DO ig=1,ngrid + transm(ig)=exp(-coef*sqrt(dup(ig))) + ENDDO + ELSE + DO ig=1,ngrid + transm(ig)=exp(-coef*dup(ig)) + ENDDO + ENDIF + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/mucorr.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/mucorr.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/mucorr.F (revision 4176) @@ -0,0 +1,107 @@ + SUBROUTINE mucorr(npts,pdeclin, plat, pmu, pfract,phaut,prad) + IMPLICIT NONE + +c======================================================================= +c +c Calcul of equivalent solar angle and and fraction of day whithout +c diurnal cycle. +c +c parmeters : +c ----------- +c +c Input : +c ------- +c npts number of points +c pdeclin solar declinaison +c plat(npts) latitude +c phaut hauteur typique de l'atmosphere +c prad rayon planetaire +c +c Output : +c -------- +c pmu(npts) equivalent cosinus of the solar angle +c pfract(npts) fractionnal day +c +c======================================================================= + +c----------------------------------------------------------------------- +c +c 0. Declarations : +c ----------------- + +c Arguments : +c ----------- + INTEGER npts + REAL plat(npts),pmu(npts), pfract(npts) + REAL phaut,prad,pdeclin +c +c Local variables : +c ----------------- + INTEGER j + REAL pi + REAL z,cz,sz,tz,phi,cphi,sphi,tphi + REAL ap,a,t,b + REAL alph + +c----------------------------------------------------------------------- + + print*,'npts,pdeclin' + print*,npts,pdeclin + pi = 4. * atan(1.0) + print*,'PI=',pi + pi=2.*asin(1.) + print*,'PI=B',pi + z = pdeclin + cz = cos (z) + sz = sin (z) + print*,'cz,sz',cz,sz + + DO 20 j = 1, npts + + phi = plat(j) + cphi = cos(phi) + if (cphi.le.1.e-9) cphi=1.e-9 + sphi = sin(phi) + tphi = sphi / cphi + b = cphi * cz + t = -tphi * sz / cz + a = 1.0 - t*t + ap = a + + IF(t.eq.0.) then + t=0.5*pi + ELSE + IF (a.lt.0.) a = 0. + t = sqrt(a) / t + IF (t.lt.0.) then + t = -atan (-t) + pi + ELSE + t = atan(t) + ENDIF + ENDIF + + pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi + pfract(j) = t / pi + IF (ap .lt.0.) then + pmu(j) = sphi * sz + pfract(j) = 1.0 + ENDIF + IF (pmu(j).le.0.0) pmu(j) = 0.0 + pmu(j) = pmu(j) / pfract(j) + IF (pmu(j).eq.0.) pfract(j) = 0. + + 20 CONTINUE + +c----------------------------------------------------------------------- +c correction de rotondite: +c ------------------------ + + alph=phaut/prad + DO 30 j=1,npts +c !!!!!! + pmu(j)=sqrt(1224.*pmu(j)*pmu(j)+1.)/35. +c $ (sqrt(alph*alph*pmu(j)*pmu(j)+2.*alph+1.)-alph*pmu(j)) +30 CONTINUE + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/mufract.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/mufract.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/mufract.F (revision 4176) @@ -0,0 +1,96 @@ + SUBROUTINE mufract(jjm,pdecli, plat, pmu, pfract) + IMPLICIT NONE +c +c======================================================================= +c +c Calcul of equivalent solar angle and and fraction of day whithout +c diurnal cycle. +c +c parmeters : +c ----------- +c +c Input : +c ------- +c jjm number of points +c pdecli solar declinaison +c plat(jjm) latitude +c +c Output : +c -------- +c pmu(jjm) equivalent cosinus of the solar angle +c pfract(jjm) fractionnal day +c +c +c inclinex is the obliquity of the planet. +c +c======================================================================= +c +c----------------------------------------------------------------------- +c +c 0. Declarations : +c ----------------- +c +c Arguments : +c ----------- + + INTEGER jjm + REAL plat(jjm),pmu(jjm), pfract(jjm) + REAL pdecli +c +c Local variables : +c ----------------- + + REAL inclinex + PARAMETER (inclinex=26.7) + INTEGER j + REAL pi + REAL z,cz,sz,tz,phi,cphi,sphi,tphi + REAL ap,a,t,b +c +c======================================================================= +c + pi = 4. * atan(1.0) + z = pdecli + cz = cos (z*pi/180.) + sz = sin (z*pi/180.) +c + DO 20 j = 1, jjm +c + phi = plat(j) + cphi = cos(phi) + if (cphi.le.1.e-9) cphi=1.e-9 + sphi = sin(phi) + tphi = sphi / cphi + b = cphi * cz + t = -tphi * sz / cz + a = 1.0 - t*t + ap = a + IF(t.eq.0.) THEN + t=0.5*pi + ELSE + IF (a.lt.0.) a = 0. + t = sqrt(a) / t + IF (t.lt.0.) THEN + t = -atan (-t) + pi + ELSE + t = atan(t) + ENDIF + ENDIF + pmu(j) = (sphi*sz*t) / pi + b*sin(t)/pi + pfract(j) = t / pi + IF (ap .lt.0.) THEN + pmu(j) = sphi * sz + pfract(j) = 1.0 + ENDIF + IF (pmu(j).le.0.0) pmu(j) = 0.0 + pmu(j) = pmu(j) / pfract(j) + IF (pmu(j).eq.0.) pfract(j) = 0. +c + 20 CONTINUE +c + RETURN + END +c +c* end of mufract +c======================================================================= +c Index: /dynamico_lmdz/simple_physics/phyparam/param/multipl.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/multipl.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/multipl.F (revision 4176) @@ -0,0 +1,17 @@ + SUBROUTINE multipl(n,x1,x2,y) + IMPLICIT NONE +c======================================================================= +c +c multiplication de deux vecteurs +c +c======================================================================= +c + INTEGER n,i + REAL x1(n),x2(n),y(n) +c + DO 10 i=1,n + y(i)=x1(i)*x2(i) +10 CONTINUE +c + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/my_25.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/my_25.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/my_25.F (revision 4176) @@ -0,0 +1,247 @@ + subroutine my_25(imax,kmax,dt,gp,zi,z,u,v,teta,cd,q2,long,km,kh) + +*********************************************************************** +******* FERMETURE MELLOR-YAMADA DE NIVEAU 2.5 (QUASI-EQUILIBRE) ******* +** q2 au interfaces entre mailles. +*********************************************************************** + + +************** DECLARATIONS ******************************************* + + + integer imax,kmax + + parameter (impmax=5) + + real kappa,khmin,kmmin,kqmin,longc + + parameter (kappa=0.4) + parameter (a1=0.92, a2=0.74, b1=16.6, b2=10.1, c1=0.08, + a e1=1.8, e2=1.33) + parameter (khmin=1.0e-5, kmmin=1.0e-5, kqmin=1.0e-5, + a q2min=0.001, q2lmin=0.001) + parameter (ghmax=0.023, ghmin=-0.28) + + real cd(imax),teta(imax,kmax),u(imax,kmax),v(imax,kmax), + a z(imax,kmax),zi(imax,kmax+1) + real kh(imax,kmax+1),km(imax,kmax+1),q2(imax,kmax+1), + a long(imax,kmax+1) + real unsdz(imax,kmax),unsdzi(imax,kmax+1) + real kq(imax,kmax), + a m2(imax,kmax+1),n2(imax,kmax+1),ri(imax,kmax+1) + real a(imax,kmax),b(imax,kmax),c(imax,kmax),f(imax,kmax), + a alph(imax,kmax) + real ksdz2inf,ksdz2sup + +c save q2,q2l + + save imp + + data imp/0/ + +*********************************************************************** + + imp=imp+1 + +************** INCREMENTS VERTICAUX *********************************** + + do 9 i=1,imax + zi(i,kmax+1)=zi(i,kmax)+2.0*(z(i,kmax)-zi(i,kmax)) + 9 continue + do 10 k=1,kmax + do 10 i=1,imax + unsdz(i,k)=1.0/(zi(i,k+1)-zi(i,k)) + 10 continue + do 11 k=2,kmax + do 11 i=1,imax + unsdzi(i,k)=1.0/(z(i,k)-z(i,k-1)) + 11 continue + do 12 i=1,imax + unsdzi(i,1)=0.5/(z(i,1)-zi(i,1)) + unsdzi(i,kmax+1)=0.5/(zi(i,kmax+1)-z(i,kmax)) + 12 continue + +*********************************************************************** + + +************** DIFFUSIVITES KH, KM et KQ ****************************** +* Ci-dessous, une premiere estimation des diffusivites turbulentes km * +* et kh est effectuee pour utilisation dans les taux de production * +* et de destruction de q2 et q2l. On calcule aussi kq. * + + do 100 k=2,kmax + do 100 i=1,imax + beta=2.0/(teta(i,k)+teta(i,k-1)) + n2(i,k)=beta*gp*unsdzi(i,k)*(teta(i,k)-teta(i,k-1)) + n2(i,k)=amax1(0.0,n2(i,k)) + du=unsdzi(i,k)*(u(i,k)-u(i,k-1)) + dv=unsdzi(i,k)*(v(i,k)-v(i,k-1)) + m2(i,k)=du*du+dv*dv + ri(i,k)=n2(i,k)/(m2(i,k)+1.0e-10) + ri(i,k)=amax1(-0.1,min(4.0,ri(i,k))) + 100 continue + + do 110 k=2,kmax + do 110 i=1,imax + vt=kappa*(zi(i,k)-zi(i,1)) + long(i,k)=vt/(1.0+vt/160.0) + if(n2(i,k).gt.0.0) then + long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) + endif + gh=amax1(ghmin, + a min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k))) + sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh* + a ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/ + a ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh)) + km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) + sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) + kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) + 110 continue + do 111 i=1,imax + us=sqrt(cd(i)*(u(i,1)*u(i,1)+v(i,1)*v(i,1))) + vt1=(b1**0.666667)*us*us + vt2=(b1**0.6666667)*kappa*kappa* + a m2(i,2)*(zi(i,2)-zi(i,1))*(zi(i,2)-zi(i,1)) +c q2(i,1)=amax1(q2min,vt1,vt2) + q2(i,1)=amax1(q2min,vt1) + long(i,1)=0.0 + long(i,kmax+1)=long(i,kmax) + sq=0.2 + kq(i,1)=amax1(kqmin,kappa*(z(i,1)-zi(i,1))*us*sq) + 111 continue + + do 120 k=2,kmax + do 120 i=1,imax + longc=0.5*(long(i,k)+long(i,k+1)) + q2c=0.5*(q2(i,k)+q2(i,k+1)) + sq=0.2 + kq(i,k)=amax1(kqmin,longc*sqrt(q2c)*sq) + 120 continue + +*********************************************************************** + + +************** CALCUL DE Q2 ******************************************* + + do 200 k=2,kmax + do 200 i=1,imax + prod=2.0*(km(i,k)*m2(i,k)+amax1(0.0,-kh(i,k)*n2(i,k))) + dest=2.0*(amax1(0.0,kh(i,k)*n2(i,k)) + a +q2(i,k)*sqrt(q2(i,k))/(b1*long(i,k))) + if(k.lt.kmax) then + ksdz2sup=unsdzi(i,k)*unsdz(i,k)*kq(i,k) + else + ksdz2sup=0.0 + endif + ksdz2inf=unsdzi(i,k)*unsdz(i,k-1)*kq(i,k-1) + b(i,k)=-ksdz2inf*dt + a(i,k)=1.0+dt*(dest/q2(i,k)+ksdz2inf+ksdz2sup) + c(i,k)=-ksdz2sup*dt + f(i,k)=q2(i,k)+dt*prod + 200 continue + do 201 i=1,imax + f(i,2)=f(i,2) + a +dt*unsdzi(i,2)*unsdz(i,1)*kq(i,1)*q2(i,1) + 201 continue + + do 210 i=1,imax + alph(i,2)=a(i,2) + 210 continue + do 211 k=3,kmax + do 211 i=1,imax + bet=b(i,k)/alph(i,k-1) + alph(i,k)=a(i,k)-bet*c(i,k-1) + f(i,k)=f(i,k)-bet*f(i,k-1) + 211 continue + + do 220 i=1,imax + q2(i,kmax)=amax1(q2min,f(i,kmax)/alph(i,kmax)) + q2(i,kmax+1)=q2(i,kmax) + 220 continue + do 221 k=kmax-1,2,-1 + do 221 i=1,imax + q2(i,k)=amax1(q2min,(f(i,k)-c(i,k)*q2(i,k+1))/alph(i,k)) + 221 continue + do 222 i=1,imax + q2(i,2)=amax1(q2(i,2),q2(i,1)) + 222 continue + +*********************************************************************** + + +************** EVALUATION FINALE DE KH ET KM ************************** + + do 400 k=2,kmax + do 400 i=1,imax + if(n2(i,k).gt.0.0) then + long(i,k)=min(0.53*sqrt(q2(i,k))/sqrt(n2(i,k)),long(i,k)) + endif + gh=amax1(ghmin, + a min(ghmax,-long(i,k)*long(i,k)*n2(i,k)/q2(i,k))) + sm=a1*(1.0-3.0*c1-6.0*a1/b1-3.0*a2*gh* + a ((b2-3.0*a2)*(1.0-6.0*a1/b1)-3.0*c1*(b2+6.0*a1)))/ + a ((1.0-3.0*a2*gh*(6.0*a1+b2))*(1.0-9.0*a1*a2*gh)) + km(i,k)=amax1(kmmin,long(i,k)*sqrt(q2(i,k))*sm) + sh=a2*(1.0-6.0*a1/b1)/(1.0-3.0*a2*gh*(6.0*a1+b2)) + kh(i,k)=amax1(khmin,long(i,k)*sqrt(q2(i,k))*sh) + 400 continue + do 401 i=1,imax + km(i,1)=kmmin + km(i,kmax+1)=km(i,kmax) + kh(i,1)=khmin + kh(i,kmax+1)=kh(i,kmax) + 401 continue + +*********************************************************************** + +c if(imp.eq.impmax) then + am=1.0/float(imax) + imp=0 + do 1000 k=kmax,1,-1 + au=0.0 + ateta=0.0 + az=0.0 + adz=0.0 + akq=0.0 + acd=0.0 + do 1001 i=1,imax + au=au+am*sqrt(u(i,k)*u(i,k)+v(i,k)*v(i,k)) + ateta=ateta+am*teta(i,k) + az=az+am*z(i,k) + adz=adz+am*(1.0/unsdz(i,k)) + akq=akq+am*kq(i,k) + acd=acd+am*cd(i) + 1001 continue +c write(*,2000) k,az,adz,au,ateta,akq,acd*1000.0 + 2000 format(2x,i3,2x,6(f9.2,2x)) + 1000 continue + + write(*,*) + write(*,*) + + do 1002 k=kmax+1,1,-1 + azi=0.0 + adzi=0.0 + aq2=0.0 + al=0.0 + akm=0.0 + akh=0.0 + am2=0.0 + al0=0.0 + do 1003 i=1,imax + azi=azi+am*zi(i,k) + adzi=adzi+am*(1.0/unsdzi(i,k)) + aq2=aq2+am*q2(i,k) + al=al+am*long(i,k) + akm=akm+am*km(i,k) + akh=akh+am*kh(i,k) + am2=am2+am*m2(i,k) +c al0=al0+am*long0d(i) + 1003 continue +c write(*,2001) k,azi,aq2,al,akm,akh,am2*1.0e5 + 2001 format(2x,i3,6(2x,f9.3)) + 1002 continue +c endif + + return + end Index: /dynamico_lmdz/simple_physics/phyparam/param/orbite.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/orbite.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/orbite.F (revision 4176) @@ -0,0 +1,52 @@ + SUBROUTINE orbite(pls,pdist_sol,pdecli) + IMPLICIT NONE + +c======================================================================= +c +c Objet: +c ------ +c +c Distance from sun and declimation as a function of the solar +c longitude Ls +c +c Interface: +c ---------- +c +c common 'comorbit.h' initialized by comorbit.f +c +c Arguments: +c ---------- +c +c Input: +c ------ +c pls Ls +c +c Output: +c ------- +c pdist_sol Distance Sun-Planet in UA +c pdecli declinaison ( en radians ) +c +c======================================================================= +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "comorbit.h" + +c arguments: +c ---------- + + REAL pday,pdist_sol,pdecli,pls + +c----------------------------------------------------------------------- + +c Distance Sun-Planet + + pdist_sol=p_elips/(1.+e_elips*cos(pls+timeperi)) + +c Solar declination + + pdecli= asin (sin(pls)*sin(obliquit*pi/180.)) + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/paramdef.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/paramdef.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/paramdef.F (revision 4176) @@ -0,0 +1,174 @@ + SUBROUTINE paramdef(ngrid,rnatur,albedo,inertie,emissiv,z0) + USE comgeomfi, ONLY : lati,sinlat,coslat + IMPLICIT NONE +c----------------------------------------------------------------------- +c declarations: +c ------------- + +#include "planete.h" + + INTEGER ngrid + REAL rnatur(ngrid) + REAL albedo(ngrid),inertie(ngrid),emissiv(ngrid) + REAL z0(ngrid) + +c local: +c ------ + + character str1*1 + INTEGER ig,ierr + real zz,z1,z2,pi + + REAL I_mer,I_ter,z0_mer,z0_ter + REAL alb_mer,alb_ter,emi_mer,emi_ter + +c----------------------------------------------------------------------- +c Caracteristiques orbitales: +c----------------------------------------------------------------------- + + pi=2.*asin(1.) + year_day=360. + periheli=173. + aphelie=173. + peri_day=0. + obliquit=24. + z0_mer=0.1 + z0_ter=10 + I_mer=30000. + I_ter=2000. + alb_ter=.112 + alb_mer=.112 + emi_ter=1. + emi_mer=1. + emin_turb=1.e-6 + lmixmin=30. + coefvis=.99 + coefir=.08 + + print*,'Do you want to change the default values for the' + print*,'physical paremetrisations (y,n)' + read(*,'(a)',iostat=ierr) str1 + do while (ierr.eq.0.and.str1.eq.'y') + print*,'Do you want to change the following parameters?' + print*,'For each parameter enter a new value or RETURN' + print*,'At the end, you will have another chance' + + print* + print*,'ASTRONOMY:' + + print*,'Length of year: ',year_day,' (sideral days)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) year_day=zz + + print*,'Distance to Sun at perihelion: ',periheli,' (Mkm)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) periheli=zz + + print*,'Distance to Sun at aphelion: ',aphelie,' (Mkm)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) aphelie=zz + + print*,'Date of the perihelion: ',peri_day,' (days)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) peri_day=zz + + print*,'Obliquity: ',obliquit,' (degrees)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) obliquit=zz + + print* + print*,'SURFACE:' + + print*,'roughness for 1. oceans and 2. solid surfaces: ', + s z0_mer,z0_ter,' (m)' + read(*,*,iostat=ierr) z1,z2 + if(ierr.eq.0) then + z0_mer=z1 + z0_ter=z2 + endif + + print*,'Thermal inertia for 1. oceans and 2. solid surfaces: ', + s I_mer,I_ter,' (SI)' + read(*,*,iostat=ierr) z1,z2 + if(ierr.eq.0) then + I_mer=z1 + I_ter=z2 + endif + + print*,'Visible albedo for 1. oceans and 2. solid surfaces: ', + s alb_mer,alb_ter + read(*,*,iostat=ierr) z1,z2 + if(ierr.eq.0) then + alb_mer=z1 + alb_ter=z2 + endif + + print*,'thermal emissivity for 1. oceans and 2. solid surf.: ', + s emi_mer,emi_ter + read(*,*,iostat=ierr) z1,z2 + if(ierr.eq.0) then + emi_mer=z1 + emi_ter=z2 + endif + + print* + print*,'PLANETARY BOUNDARY LAYER:' + + print*,'Minimum turbulent kinetic energie: ',emin_turb, + s ' (m2s-2)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) emin_turb=zz + + print*,'Mixing length: ',lmixmin,' (m)' + read(*,*,iostat=ierr) zz + if(ierr.eq.0) lmixmin=zz + + print* + print*,'RADIATION:' + + print*,'Fraction of the incoming solar radiation ', + s 'reaching the surface directly: ',coefvis + read(*,*,iostat=ierr) zz + if(ierr.eq.0) coefvis=zz + + print*,'Fraction of the surface thermal emission ', + s 'reaching space directly: ',coefir + read(*,*,iostat=ierr) zz + if(ierr.eq.0) coefir=zz + + print*,'Do you want to change something once again?' + read(*,'(a)',iostat=ierr) str1 + enddo + +c----------------------------------------------------------------------- + + DO ig=1,ngrid +c inertie(ig)=(1.-rnatur(ig))*I_mer+rnatur(ig)*I_ter +c albedo(ig)=(1.-rnatur(ig))*alb_mer+rnatur(ig)*alb_ter + if (rnatur(ig).gt.0.5) then + if (lati(ig).lt.-pi/3.) then +c antartique + inertie(ig)=100000. + albedo(ig)=0.6 + else +c continents normaux + inertie(ig)=2000. + albedo(ig)=0.112 + endif + else +c oceans + inertie(ig)=100000. + albedo(ig)=0.112 + endif + + z0(ig)=(1.-rnatur(ig))*z0_mer+rnatur(ig)*z0_ter + emissiv(ig)=(1.-rnatur(ig))*emi_mer+rnatur(ig)*emi_ter + if(rnatur(ig).ge..5.and.coslat(ig).lt..5.and.sinlat(ig).lt.0.) + . albedo(ig)=.5 + ENDDO + + RETURN + +9000 FORMAT(3(/,a72)) +9001 FORMAT(/,a72,/,a8) + END Index: /dynamico_lmdz/simple_physics/phyparam/param/phyaqua.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/phyaqua.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/phyaqua.F (revision 4176) @@ -0,0 +1,26 @@ +! Routines complementaires pour la physique planetaire. + + + subroutine iniaqua(nlon,latfi,lonfi,iflag_phys) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Creation d'un etat initial et de conditions aux limites +! (resp startphy.nc et limit.nc) pour des configurations idealisees +! du modele LMDZ dans sa version terrestre. +! iflag_phys est un parametre qui controle +! iflag_phys = N +! de 100 a 199 : aqua planetes avec SST forcees +! N-100 determine le type de SSTs +! de 200 a 299 : terra planetes avec Ts calcule +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + integer nlon,iflag_phys +cIM ajout latfi, lonfi + REAL, DIMENSION (nlon) :: lonfi, latfi + + + return + end + Index: /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/phyparam.F (revision 4176) @@ -0,0 +1,661 @@ + SUBROUTINE phyparam(ngrid,nlayer,nq, + s firstcall,lastcall, + s rjourvrai,gmtime,ptimestep, + s pplev,pplay,pphi,pphis,presnivs, + s pu,pv,pt,pq, + s pw, + s pdu,pdv,pdt,pdq,pdpsrf) + + USE comsaison + USE dimphy + USE comgeomfi +c + IMPLICIT NONE +c======================================================================= +c +c subject: +c -------- +c +c Organisation of the physical parametrisations of the LMD +c 20 parameters GCM for planetary atmospheres. +c It includes: +c raditive transfer (long and shortwave) for CO2 and dust. +c vertical turbulent mixing +c convective adjsutment +c +c author: Frederic Hourdin 15 / 10 /93 +c ------- +c +c arguments: +c ---------- +c +c input: +c ------ +c +c ngrid Size of the horizontal grid. +c All internal loops are performed on that grid. +c nlayer Number of vertical layers. +c nq Number of advected fields +c firstcall True at the first call +c lastcall True at the last call +c rjourvrai Number of days counted from the North. Spring +c equinoxe. +c gmtime hour (s) +c ptimestep timestep (s) +c pplay(ngrid,nlayer+1) Pressure at the middle of the layers (Pa) +c pplev(ngrid,nlayer+1) intermediate pressure levels (pa) +c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2) +c pu(ngrid,nlayer) u component of the wind (ms-1) +c pv(ngrid,nlayer) v component of the wind (ms-1) +c pt(ngrid,nlayer) Temperature (K) +c pq(ngrid,nlayer,nq) Advected fields +c pudyn(ngrid,nlayer) \ +c pvdyn(ngrid,nlayer) \ Dynamical temporal derivative for the +c ptdyn(ngrid,nlayer) / corresponding variables +c pqdyn(ngrid,nlayer,nq) / +c pw(ngrid,?) vertical velocity +c +c output: +c ------- +c +c pdu(ngrid,nlayermx) \ +c pdv(ngrid,nlayermx) \ Temporal derivative of the corresponding +c pdt(ngrid,nlayermx) / variables due to physical processes. +c pdq(ngrid,nlayermx) / +c pdpsrf(ngrid) / +c +c======================================================================= +c +c----------------------------------------------------------------------- +c +c 0. Declarations : +c ------------------ + +#include "dimensions.h" +#include "description.h" +#include "callkeys.h" +#include "comcstfi.h" +#include "planete.h" +#include "surface.h" + +c Arguments : +c ----------- + +c inputs: +c ------- + INTEGER ngrid,nlayer,nq + + REAL ptimestep + real zdtime + REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) + REAL pphi(ngrid,nlayer) + REAL pphis(ngrid) + REAL pu(ngrid,nlayer),pv(ngrid,nlayer) + REAL pt(ngrid,nlayer),pq(ngrid,nlayer,nq) + REAL pdu(ngrid,nlayer),pdv(ngrid,nlayer) + +c dynamial tendencies + REAL pdtdyn(ngrid,nlayer),pdqdyn(ngrid,nlayer,nq) + REAL pdudyn(ngrid,nlayer),pdvdyn(ngrid,nlayer) + REAL pw(ngrid,nlayer) + +c Time + real rjourvrai + REAL gmtime + +c outputs: +c -------- + +c physical tendencies + REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq) + REAL pdpsrf(ngrid) + LOGICAL firstcall,lastcall + +c Local variables : +c ----------------- + + INTEGER j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq + INTEGER*4 day_ini +c + REAL,DIMENSION(ngrid) :: mu0,fract + REAL zday + REAL zh(ngrid,nlayer),z1,z2 + REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer) + REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer) + REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid) + REAL zflubid(ngrid),zpmer(ngrid) + REAL zplanck(ngrid),zpopsk(ngrid,nlayer) + REAL zdum1(ngrid,nlayer) + REAL zdum2(ngrid,nlayer) + REAL zdum3(ngrid,nlayer) + REAL ztim1,ztim2,ztim3 + REAL zls,zinsol + REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer) + REAL zfluxsw(ngrid),zfluxlw(ngrid) + character*2 str2 + REAL factq(nq),tauq(nq) + character*3 nomq + +c Local saved variables: +c ---------------------- + + INTEGER, SAVE :: icount + real, SAVE :: zday_last +!$OMP THREADPRIVATE( icount,zday_last) + + REAL zps_av + + real,allocatable,SAVE :: tsurf(:),tsoil(:,:),rnatur(:) + real,allocatable,SAVE :: capcal(:),fluxgrd(:) + real,allocatable,SAVE :: dtrad(:,:),fluxrad(:) + real,allocatable,SAVE :: q2(:,:),q2l(:,:) + real,allocatable,SAVE :: albedo(:),emissiv(:),z0(:),inertie(:) + real,SAVE :: solarcst=1370. + real,SAVE :: stephan=5.67e-08 + +!$OMP THREADPRIVATE(tsurf,tsoil,rnatur) +!$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad) +!$OMP THREADPRIVATE( q2,q2l) +!$OMP THREADPRIVATE( albedo,emissiv,solarcst,z0,inertie) +!$OMP THREADPRIVATE( stephan) + + + EXTERNAL vdif,convadj + EXTERNAL orbite,mucorr + EXTERNAL ismin,ismax + + + INTEGER longcles + PARAMETER ( longcles = 20 ) + REAL clesphy0( longcles ) + REAL presnivs(nlayer) + + + + print*,'OK DANS PHYPARAM' + +c----------------------------------------------------------------------- +c 1. Initialisations : +c -------------------- + +! call initial0(ngrid*nlayermx*nqmx,pdq) + nlevel=nlayer+1 + +! print*,'OK ',nlevel + + igout=ngrid/2+1 +! print*,'OK PHYPARAM ',ngrid,igout + + + zday=rjourvrai+gmtime + +! print*,'OK PHYPARAM 0A nq ',nq + tauq(1)=1800. + tauq(2)=10800. + tauq(3)=86400. + tauq(4)=864000. +! print*,'OK PHYPARAM 0 B nq ',nq + factq(1:4)=(1.-exp(-ptimestep/tauq(1:4)))/ptimestep + +! print*,'OK PHYPARAM 0 ' + print*,'nq ',nq + print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid) + print*,'nlayer',nlayer + print*,'size pdq ',ngrid*nlayer*4,ngrid*nlayer*nq, + s size(pdq),size(lati),size(pq),size(factq) + do iq=1,4 + do l=1,nlayer + pdq(1:ngrid,l,iq)= + & (1.+sin(lati(1:ngrid))-pq(1:ngrid,l,iq))*factq(iq) + enddo + enddo + + IF(firstcall) THEN + + print*,'AKk',ngrid,nsoilmx + allocate(tsurf(ngrid)) + print*,'AKa' + allocate (tsoil(ngrid,nsoilmx)) + print*,'AKb' + allocate (rnatur(ngrid)) + print*,'AK2' + allocate(capcal(ngrid),fluxgrd(ngrid)) + print*,'AK3' + allocate(dtrad(ngrid,nlayer),fluxrad(ngrid)) + print*,'AK4' + allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1)) + print*,'AK5' + allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid)) + print*,'AK6' + + + do l=1,nlayer + pdq(:,l,5)=1.+sin(lati(:))/ptimestep + enddo + PRINT*,'FIRSTCALL ' + +! zday_last=rjourvrai + zday_last=zday-ptimestep/unjours +c CALL readfi(ngrid,nlayer,nsoilmx,ldrs, +c . day_ini,gmtime,albedo,inertie,emissiv,z0,rnatur, +c . q2,q2l,tsurf,tsoil) + rnatur=1. + emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter + inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter + albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter + z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter + q2=1.e-10 + q2l=1.e-10 + tsurf=300. + tsoil=300. + print*,tsoil(ngrid/2+1,nsoilmx/2+2) + print*,'TS ',tsurf(igout),tsoil(igout,5) + CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) + + if (.not.callrad) then + do ig=1,ngrid + fluxrad(ig)=0. + enddo + endif + +! print*,'OK PHYPARAM 1 ' + IF(callsoil) THEN + CALL soil(ngrid,nsoilmx,firstcall,inertie, + s ptimestep,tsurf,tsoil,capcal,fluxgrd) + ELSE + PRINT*,'WARNING!!! Thermal conduction in the soil + s turned off' + DO ig=1,ngrid + capcal(ig)=1.e5 + fluxgrd(ig)=0. + ENDDO + ENDIF +c CALL inifrict(ptimestep) + icount=0 + + PRINT*,'FIRSTCALL B ' + print*,'INIIO avant iophys_ini ' + call iophys_ini('phys.nc ',presnivs) + + ENDIF + icount=icount+1 + + PRINT*,'FIRSTCALL AP ' + IF (ngrid.NE.ngridmax) THEN + PRINT*,'STOP in inifis' + PRINT*,'Probleme de dimenesions :' + PRINT*,'ngrid = ',ngrid + PRINT*,'ngridmax = ',ngridmax + STOP + ENDIF + + DO l=1,nlayer + DO ig=1,ngrid + pdv(ig,l)=0. + pdu(ig,l)=0. + pdt(ig,l)=0. + ENDDO + ENDDO + + DO ig=1,ngrid + pdpsrf(ig)=0. + zflubid(ig)=0. + zdtsrf(ig)=0. + ENDDO + +c----------------------------------------------------------------------- +c calcul du geopotentiel aux niveaux intercouches +c ponderation des altitudes au niveau des couches en dp/p + + DO l=1,nlayer + DO ig=1,ngrid + zzlay(ig,l)=pphi(ig,l)/g + ENDDO + ENDDO + DO ig=1,ngrid + zzlev(ig,1)=0. + ENDDO + DO l=2,nlayer + DO ig=1,ngrid + z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) + z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) + zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c Transformation de la temperature en temperature potentielle + DO l=1,nlayer + DO ig=1,ngrid + zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp + ENDDO + ENDDO + DO l=1,nlayer + DO ig=1,ngrid + zh(ig,l)=pt(ig,l)/zpopsk(ig,l) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 2. Calcul of the radiative tendencies : +c --------------------------------------- + +! print*,'callrad0' + IF(callrad) THEN +! print*,'callrad' + +c WARNING !!! on calcule le ray a chaque appel +c IF( MOD(icount,iradia).EQ.0) THEN + + CALL solarlong(zday,zls) + CALL orbite(zls,dist_sol,declin) +c declin=0. +! print*,'ATTENTIOn : pdeclin = 0',' L_s=',zls + print*,'diurnal=',diurnal + IF(diurnal) THEN + if ( 1.eq.1 ) then + ztim1=SIN(declin) + ztim2=COS(declin)*COS(2.*pi*(zday-.5)) + ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) +c call dump2d(iim,jjm-1,lati(2),'LATI 0 ') +c call dump2d(iim,jjm-1,long(2),'LONG 0 ') +c call dump2d(iim,jjm-1,sinlon(2),'sinlon0 ') +c call dump2d(iim,jjm-1,coslon(2),'coslon0 ') +c call dump2d(iim,jjm-1,sinlat(2),'sinlat ') +c call dump2d(iim,jjm-1,coslat(2),'coslat ') + CALL solang(ngrid,sinlon,coslon,sinlat,coslat, + s ztim1,ztim2,ztim3, + s mu0,fract) + else + zdtime=ptimestep*float(iradia) + CALL zenang(zls,gmtime,zdtime,lati,long,mu0,fract) + print*,'ZENANG ' + endif + + IF(lverbose) THEN + PRINT*,'day, declin, sinlon,coslon,sinlat,coslat' + PRINT*,zday, declin, + s sinlon(igout),coslon(igout), + s sinlat(igout),coslat(igout) + ENDIF + ELSE + print*,'declin,ngrid,rad',declin,ngrid,rad + +c call dump2d(iim,jjm-1,lati(2),'LATI ') + CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad) + ENDIF +c call dump2d(iim,jjm-1,fract(2),'FRACT A ') +c call dump2d(iim,jjm-1,mu0(2),'MU0 A ') + + +c 2.2 Calcul of the radiative tendencies and fluxes: +c -------------------------------------------------- + +c 2.1.2 levels + + zinsol=solarcst/(dist_sol*dist_sol) + print*,iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer + print*,'iim,jjm,llm,ngrid,nlayer,ngridmax,nlayer' +c call dump2d(iim,jjm-1,albedo(2),'ALBEDO ') +c call dump2d(iim,jjm-1,mu0(2),'MU0 ') +c call dump2d(iim,jjm-1,fract(2),'FRACT ') +c call dump2d(iim,jjm-1,lati(2),'LATI ') + zps_av=1.e5 + if (firstcall) then + print*,'WARNING ps_rad impose' + endif + CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, + $ pplev,zps_av, + $ mu0,fract,zinsol, + $ zfluxsw,zdtsw, + $ lverbose) +c call dump2d(iim,jjm-1,zfluxsw(2),'SWS 1 ') +c stop + + CALL lw(ngrid,nlayer,coefir,emissiv, + $ pplev,zps_av,tsurf,pt, + $ zfluxlw,zdtlw, + $ lverbose) + + +c 2.4 total flux and tendencies: +c ------------------------------ + +c 2.4.1 fluxes + + DO ig=1,ngrid + fluxrad(ig)=emissiv(ig)*zfluxlw(ig) + $ +zfluxsw(ig)*(1.-albedo(ig)) + zplanck(ig)=tsurf(ig)*tsurf(ig) + zplanck(ig)=emissiv(ig)* + $ stephan*zplanck(ig)*zplanck(ig) + fluxrad(ig)=fluxrad(ig)-zplanck(ig) + ENDDO + +c 2.4.2 temperature tendencies + + DO l=1,nlayer + DO ig=1,ngrid + dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) + ENDDO + ENDDO + +c ENDIF + +c 2.5 Transformation of the radiative tendencies: +c ----------------------------------------------- + + DO l=1,nlayer + DO ig=1,ngrid + pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) + ENDDO + ENDDO + + IF(lverbose) THEN + PRINT*,'Diagnotique for the radaition' + PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck' + PRINT*,albedo(igout),emissiv(igout),mu0(igout), + s fract(igout), + s fluxrad(igout),zplanck(igout) + PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)' + PRINT*,'unjours',unjours + DO l=1,nlayer + WRITE(*,'(3f15.5,2e15.2)') pt(igout,l), + s pplay(igout,l),pplev(igout,l), + s zdtsw(igout,l),zdtlw(igout,l) + ENDDO + ENDIF + + + ENDIF + +c----------------------------------------------------------------------- +c 3. Vertical diffusion (turbulent mixing): +c ----------------------------------------- +c + IF(calldifv) THEN + + DO ig=1,ngrid + zflubid(ig)=fluxrad(ig)+fluxgrd(ig) + ENDDO + + CALL zerophys(ngrid*nlayer,zdum1) + CALL zerophys(ngrid*nlayer,zdum2) +c CALL zerophys(ngrid*nlayer,zdum3) + do l=1,nlayer + do ig=1,ngrid + zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l) + enddo + enddo + + CALL vdif(ngrid,nlayer,zday, + $ ptimestep,capcal,z0, + $ pplay,pplev,zzlay,zzlev, + $ pu,pv,zh,tsurf,emissiv, + $ zdum1,zdum2,zdum3,zflubid, + $ zdufr,zdvfr,zdhfr,zdtsrfr,q2,q2l, + $ lverbose) +c igout=ngrid/2+1 +c PRINT*,'zdufr zdvfr zdhfr' +c DO l=1,nlayer +c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l) +c ENDDO +c CALL difv (ngrid,nlayer, +c $ area,lati,capcal, +c $ pplev,pphi, +c $ pu,pv,zh,tsurf,emissiv, +c $ zdum1,zdum2,zdum3,zflubid, +c $ zdufr,zdvfr,zdhfr,zdtsrf, +c $ .true.) +c PRINT*,'zdufr zdvfr zdhfr' +c DO l=1,nlayer +c PRINT*,zdufr(igout,1),zdvfr(igout,l),zdhfr(igout,l) +c ENDDO +c STOP + + DO l=1,nlayer + DO ig=1,ngrid + pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) + pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) + pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) + ENDDO + ENDDO + + DO ig=1,ngrid + zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig) + ENDDO + + ELSE + DO ig=1,ngrid + zdtsrf(ig)=zdtsrf(ig)+ + s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) + ENDDO + ENDIF +c +c----------------------------------------------------------------------- +c 4. Dry convective adjustment: +c ----------------------------- + + IF(calladj) THEN + + DO l=1,nlayer + DO ig=1,ngrid + zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l) + ENDDO + ENDDO + CALL zerophys(ngrid*nlayer,zdufr) + CALL zerophys(ngrid*nlayer,zdvfr) + CALL zerophys(ngrid*nlayer,zdhfr) + CALL convadj(ngrid,nlayer,ptimestep, + $ pplay,pplev,zpopsk, + $ pu,pv,zh, + $ pdu,pdv,zdum1, + $ zdufr,zdvfr,zdhfr) + + DO l=1,nlayer + DO ig=1,ngrid + pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) + pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) + pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) + ENDDO + ENDDO + + ENDIF + +c----------------------------------------------------------------------- +c On incremente les tendances physiques de la temperature du sol: +c --------------------------------------------------------------- + + DO ig=1,ngrid + tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) + ENDDO + + WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1) + +c----------------------------------------------------------------------- +c soil temperatures: +c -------------------- + + IF (callsoil) THEN + CALL soil(ngrid,nsoilmx,.false.,inertie, + s ptimestep,tsurf,tsoil,capcal,fluxgrd) + IF(lverbose) THEN + PRINT*,'Surface Heat capacity,conduction Flux, Ts, + s dTs, dt' + PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), + s zdtsrf(igout),ptimestep + ENDIF + ENDIF + +c----------------------------------------------------------------------- +c sorties: +c -------- + +c call dump2d(iim,jjm-1,zfluxsw(2),'SWS 2 ') + print*,'zday, zday_last ',zday,zday_last,icount + if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then + print*,'zday, zday_last SORTIE ',zday,zday_last + zday_last=zday +c Ecriture/extension de la coordonnee temps + + do ig=1,ngridmax + zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*285.)) + enddo + + call iophys_ecrit('u',llm,'Vent zonal moy','m/s',pu) + call iophys_ecrit('v',llm,'Vent meridien moy','m/s',pv) + call iophys_ecrit('temp',llm,'Temperature','K',pt) + call iophys_ecrit('geop',llm,'Geopotential','m2/s2',pphi) + call iophys_ecrit('plev',llm,'plev','Pa',pplev(:,1:nlayer)) + + call iophys_ecrit('du',llm,'du',' ',pdu) + call iophys_ecrit('dv',llm,'du',' ',pdv) + call iophys_ecrit('dt',llm,'du',' ',pdt) + call iophys_ecrit('dtsw',llm,'dtsw',' ',zdtsw) + call iophys_ecrit('dtlw',llm,'dtlw',' ',zdtlw) + + do iq=1,nq + nomq="tr." + write(nomq(2:3),'(i1.1)') iq + call iophys_ecrit(nomq,llm,nomq,' ',pq(:,:,iq)) + enddo + + call iophys_ecrit('ts',1,'Surface temper','K',tsurf) + call iophys_ecrit('coslon',1,'coslon',' ',coslon) + call iophys_ecrit('sinlon',1,'sinlon',' ',sinlon) + call iophys_ecrit('coslat',1,'coslat',' ',coslat) + call iophys_ecrit('sinlat',1,'sinlat',' ',sinlat) + call iophys_ecrit('mu0',1,'mu0',' ',mu0) + call iophys_ecrit('alb',1,'alb',' ',albedo) + call iophys_ecrit('fract',1,'fract',' ',fract) + call iophys_ecrit('ps',1,'Surface pressure','Pa',pplev) + call iophys_ecrit('slp',1,'Sea level pressure','Pa',zpmer) + call iophys_ecrit('swsurf',1,'SW surf','Pa',zfluxsw) + call iophys_ecrit('lwsurf',1,'LW surf','Pa',zfluxlw) + + endif + +c----------------------------------------------------------------------- + IF(lastcall) THEN + call iotd_fin + PRINT*,'Ecriture du fichier de reinitialiastion de la physique' +! if (ierr.ne.0) then +! write(6,*)' Pb d''ouverture du fichier restart' +! write(6,*)' ierr = ', ierr +! call exit(1) +! endif + write(75) tsurf,tsoil +c s (tsurf(1),ig=1,iim+1), +c s ( (tsurf(ig),ig=(j-2)*iim+2,(j-1)*iim+1), +c s tsurf((j-2)*iim+2) ,j=2,jjm), +c s (tsurf(ngridmax),ig=1,iim+1), +c s ( (tsoil(1,l),ig=1,iim+1), +c s ( (tsoil(ig,l),ig=(j-2)*iim+2,(j-1)*iim+1), +c s tsoil((j-2)*iim+2,l) ,ig=2,jjm), +c s (tsoil(ngridmax,l),ig=1,iim+1) +c s ,l=1,nsoilmx) + ENDIF + + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/planete.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/planete.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/planete.h (revision 4176) @@ -0,0 +1,14 @@ +c----------------------------------------------------------------------- +c INCLUDE planet.h + + COMMON/planet/aphelie,periheli,year_day,peri_day, + $ obliquit, + $ lmixmin,emin_turb, + $ coefvis,coefir + + REAL aphelie,periheli,year_day,peri_day, + $ obliquit, + $ lmixmin,emin_turb, + $ coefvis,coefir + +c----------------------------------------------------------------------- Index: /dynamico_lmdz/simple_physics/phyparam/param/radia1d.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/radia1d.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/radia1d.F (revision 4176) @@ -0,0 +1,250 @@ + PROGRAM radia1d + USE dimphy + USE comsaison + USE comgeomfi + IMPLICIT NONE + +#ifdef TOTO + +#include "dimensions.h" +#include "surfdat.h" +#include "callkeys.h" +#include "comcstfi.h" +#include "planete.h" + +c Arguments : +c ----------- + + INTEGER ngrid,nlayer,nq + + REAL plev(nlayermx+1),play(nlayermx) + REAL temp(nlayermx) + REAL psrf + +c dynamial tendencies + + INTEGER l,ierr,aslun,nlevel,iaer +c + REAL mu0,fract + REAL day_ini,time,longitude,latitude + REAL zlay(nlayermx) + REAL ztlev(nlayermx+1) + REAL zplanck + REAL zrad(nlayermx),zradc(nlayermx) + REAL zdum1(nlayermx) + REAL zdum2(nlayermx) + REAL zdum3(nlayermx) + REAL zdum4(nlayermx) + REAL zdum5(nlayermx) + REAL stephan + REAL ztim1,ztim2,ztim3 + REAL zco2,zp + REAL ls,zmax,tau,tau_tot + REAL dtlw(nlayermx),dtsw(nlayermx) + REAL dtlwcl(nlayermx),dtswcl(nlayermx) + REAL zflux(6) + + REAL aerosol(nlayermx,5),cst_aer,pview + REAL tsurf,tsoil(nsoilmx) + REAL co2ice,albedo(2),emis + REAL dtrad(nlayermx),fluxrad + + DATA stephan/5.67e-08/ + DATA psrf,zmax/775.,9./ + DATA ls,time,latitude,longitude/0.,12.,0.,0./ + DATA diurnal/.true./ + DATA tau/.5/ + +c WARNING declin and dist_sol are prescribed instead of Ls + DATA declin,dist_sol/-24.8,1.4/ + +c----------------------------------------------------------------------- +c 1. Initialisations : +c -------------------- + + PRINT*,'tau?' + READ(5,*) tau + PRINT*,'time 0 to 24 h' + READ(5,*) time + PRINT*,'latitude (degrees)' + READ(5,*) latitude + + pi=2.*asin(1.) + ls=ls*pi/180. + time=time/24. + latitude=latitude*pi/180. + longitude=longitude*pi/180. + declin=declin*pi/180. + + ngrid=1 + nlayer=nlayermx + nlevel=nlayer+1 + + rad=3397200. + omeg=4.*asin(1.)/(88775.) + g=3.72 + mugaz=43.49 + rcp=.256793 + unjours=88775. + r = 8.314511*1000./mugaz + cpp = r/rcp + PRINT*,'Cp = ',cpp + PRINT*,'R = ',r + + CALL paramdef + + tsurf=200. + DO l=1,nlayer + zlay(l)=zmax*(l-.5)/nlayer + plev(l)=psrf*exp(-zmax*(l-1.)/nlayer) + play(l)=psrf*exp(-zlay(l)) + ENDDO + plev(nlevel)=0. + + DO l=1,nlayer + temp(l)=200. + ENDDO + + DO l=1,nlayer + dtrad(l)=0. + ENDDO + fluxrad=0. + + albedo(1)=.24 + albedo(2)=.24 + emis=1. + CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) + +c a total otempical detemph of .2 for Ps=700Pa + cst_aer=tau/700. + pview=1.66 + CALL sucst(66,19900101,8000000,g,rad,mugaz,rcp,unjours) + CALL surad(6) + + IF (ngrid.NE.ngridmax) THEN + PRINT*,'STOP in inifis' + PRINT*,'Probleme de dimenesions :' + PRINT*,'ngrid = ',ngrid + PRINT*,'ngridmax = ',ngridmax + STOP + ENDIF + +c----------------------------------------------------------------------- +c 2. Calcul of the radiative tendencies : +c --------------------------------------- + + +c CALL orbite(ls,dist_sol,declin) + IF(diurnal) THEN + ztim1=SIN(declin) + ztim2=COS(declin)*COS(2.*pi*(time-.5)) + ztim3=-COS(declin)*SIN(2.*pi*(time-.5)) + CALL solang(ngrid,sin(longitude),cos(longitude), + s sin(latitude),cos(latitude), + s ztim1,ztim2,ztim3, + s mu0,fract) + PRINT*,'time, declin, sinlon,coslon,sinlat,coslat' + PRINT*,time, declin,sin(longitude),cos(longitude), + s sin(latitude),cos(latitude) + ELSE + CALL mucorr(ngrid,declin, lati, mu0, fract,10000.,rad) + ENDIF + +c 2.2 Calcul of the radiative tendencies and fluxes: +c -------------------------------------------------- + +c 2.1.1 layers + tau_tot=0. + DO l=1,nlayer + zp=plev(1)/play(l) + aerosol(l,1)=amax1(cst_aer* + s (plev(l)-plev(l+1)) + s *exp(.03*(1.-zp)),1.e-33) + aerosol(l,1)=amax1( + s tau*(plev(l)-plev(l+1))/plev(1) + s *exp(.007*(1.-zp)),1.e-33) + tau_tot=tau_tot+aerosol(l,1) + ENDDO + PRINT*,'tau total',tau_tot + DO l=1,nlayermx + aerosol(l,1)=aerosol(l,1)*tau/tau_tot + ENDDO + +c 2.1.2 levels + +c Extrapolation for the air temperature above the surface + ztlev(1)=temp(1)+ + s (plev(1)-play(1))* + s (temp(1)-temp(2))/(play(1)-play(2)) + + DO l=2,nlevel-1 + ztlev(l)=0.5*(temp(l-1)+temp(l)) + ENDDO + + ztlev(nlevel)=temp(nlayer) + PRINT*,'temp' + + DO l=1,nlayer + zdum1(l)=0. + zdum2(l)=0. + zdum3(l)=0. + zdum4(l)=0. + zdum5(l)=0. + DO iaer=2,5 + aerosol(l,iaer)=0. + ENDDO + ENDDO + +c 2.3 Calcul of the radiative tendencies and fluxes: +c -------------------------------------------------- + + zco2=.95 + CALL radite(ngrid,nlayer,0,6,1,1, + $ aerosol,albedo,zco2,zdum4,emis, + $ mu0,zdum5, + $ plev,play,zdum1,zdum2,zdum3,ztlev,tsurf,temp,pview, + $ dtlw,dtsw,dtlwcl,dtswcl,zflux,zrad,zradc, + $ fract,dist_sol) + PRINT*,'radite' + +c 2.4 total flux and tendencies: +c ------------------------------ + +c 2.4.1 fluxes + + fluxrad=emis*zflux(4) + $ +zflux(5)*(1.-albedo(1)) + $ +zflux(6)*(1.-albedo(2)) + zplanck=tsurf*tsurf + zplanck=emis* + $ stephan*zplanck*zplanck + fluxrad=fluxrad-zplanck + +c 2.4.2 temperature tendencies + + DO l=1,nlayer + dtrad(l)=(dtsw(l)+dtlw(l))/unjours + dtsw(l)=cpp*dtsw(l)/unjours + dtlw(l)=dtlw(l) + ENDDO + + +c 2.5 Transformation of the radiative tendencies: +c ----------------------------------------------- + + PRINT*,'Diagnotique for the radaition' + PRINT*,'albedo, emissiv, mu0,fract,pview,Frad,Planck' + PRINT*,albedo(1),emis,mu0,fract,pview,fluxrad,zplanck + PRINT*,'Tlay Tlev Play Plev z aerosol dT/dt SW dT/dt LW (K/day)' + PRINT*,'unjours',unjours + DO l=1,nlayer + WRITE(*,'(2f7.1,2e11.3,f6.1,3e14.5)') + s temp(l),ztlev(l),play(l),plev(l),zlay(l), + s g*aerosol(l,1)/(plev(l)-plev(l+1)),dtsw(l),dtlw(l) + WRITE(55,'(2e15.5)') -dtlw(l),.001*zlay(l)*r*200./g + WRITE(56,'(2e15.5)') dtsw(l),zlay(l) + ENDDO + +#endif + + END Index: /dynamico_lmdz/simple_physics/phyparam/param/scatter.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/scatter.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/scatter.F (revision 4176) @@ -0,0 +1,21 @@ + subroutine monscatter(n,a,index,b) +C + implicit none + integer N,INDEX(n),I + real A(n),B(n) +c +c + DO 100 I=1,N + A(INDEX(I))=B(I) +100 CONTINUE +c + return + end +c +c +c peut etre vec_scatter(v,rindices,count,r) en f77hp +c equivalent code: +c do 99 i=1,count +c99 r(rindices(i))=v(i) +c voir p. 11-13 +c Index: /dynamico_lmdz/simple_physics/phyparam/param/soil.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/soil.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/soil.F (revision 4176) @@ -0,0 +1,193 @@ + SUBROUTINE soil(ngrid,nsoil,firstcall,ptherm_i, + s ptimestep,ptsrf,ptsoil, + s pcapcal,pfluxgrd) + IMPLICIT NONE + +c======================================================================= +c +c Auteur: Frederic Hourdin 30/01/92 +c ------- +c +c objet: computation of : the soil temperature evolution +c ------ the surfacic heat capacity "Capcal" +c the surface conduction flux pcapcal +c +c +c Method: implicit time integration +c ------- +c Consecutive ground temperatures are related by: +c T(k+1) = C(k) + D(k)*T(k) (1) +c the coefficients C and D are computed at the t-dt time-step. +c Routine structure: +c 1)new temperatures are computed using (1) +c 2)C and D coefficients are computed from the new temperature +c profile for the t+dt time-step +c 3)the coefficients A and B are computed where the diffusive +c fluxes at the t+dt time-step is given by +c Fdiff = A + B Ts(t+dt) +c or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt +c with F0 = A + B (Ts(t)) +c Capcal = B*dt +c +c Interface: +c ---------- +c +c Arguments: +c ---------- +c ngird number of grid-points +c ptimestep physical timestep (s) +c pto(ngrid,nsoil) temperature at time-step t (K) +c ptn(ngrid,nsoil) temperature at time step t+dt (K) +c pcapcal(ngrid) specific heat (W*m-2*s*K-1) +c pfluxgrd(ngrid) surface diffusive flux from ground (Wm-2) +c +c======================================================================= +c declarations: +c ------------- + + +c----------------------------------------------------------------------- +c arguments +c --------- + + INTEGER ngrid,nsoil + REAL ptimestep + REAL ptsrf(ngrid),ptsoil(ngrid,nsoil),ptherm_i(ngrid) + REAL pcapcal(ngrid),pfluxgrd(ngrid) + LOGICAL firstcall + + +c----------------------------------------------------------------------- +c local arrays +c ------------ + + INTEGER ig,jk + REAL za(ngrid),zb(ngrid) + REAL zdz2(nsoil),z1(ngrid) + REAL min_period,dalph_soil + +c local saved variables: +c ---------------------- + REAL,SAVE :: lambda + REAL,ALLOCATABLE,SAVE :: dz1(:),dz2(:),zc(:,:),zd(:,:) +!$OMP THREADPRIVATE(dz1,dz2,zc,zd,lambda) + +c----------------------------------------------------------------------- +c Depthts: +c -------- + + REAL fz,rk,fz1,rk1,rk2 + fz(rk)=fz1*(dalph_soil**rk-1.)/(dalph_soil-1.) + + print*,'firstcall soil ',firstcall + IF (firstcall) THEN + +c----------------------------------------------------------------------- +c ground levels +c grnd=z/l where l is the skin depth of the diurnal cycle: +c -------------------------------------------------------- + + print*,'nsoil,ngrid,firstcall=',nsoil,ngrid,firstcall + ALLOCATE(dz1(nsoil),dz2(nsoil)) + ALLOCATE(zc(ngrid,nsoil),zd(ngrid,nsoil)) + + min_period=20000. + dalph_soil=2. + + OPEN(99,file='soil.def',status='old',form='formatted',err=9999) + READ(99,*) min_period + READ(99,*) dalph_soil + PRINT*,'Discretization for the soil model' + PRINT*,'First level e-folding depth',min_period, + s ' dalph',dalph_soil + CLOSE(99) +9999 CONTINUE + +c la premiere couche represente un dixieme de cycle diurne + fz1=sqrt(min_period/3.14) + + DO jk=1,nsoil + rk1=jk + rk2=jk-1 + dz2(jk)=fz(rk1)-fz(rk2) + ENDDO + DO jk=1,nsoil-1 + rk1=jk+.5 + rk2=jk-.5 + dz1(jk)=1./(fz(rk1)-fz(rk2)) + ENDDO + lambda=fz(.5)*dz1(1) + PRINT*,'full layers, intermediate layers (secoonds)' + DO jk=1,nsoil + rk=jk + rk1=jk+.5 + rk2=jk-.5 + PRINT*,fz(rk1)*fz(rk2)*3.14, + s fz(rk)*fz(rk)*3.14 + ENDDO + +c Initialisations: +c ---------------- + + ELSE +c----------------------------------------------------------------------- +c Computation of the soil temperatures using the Cgrd and Dgrd +c coefficient computed at the previous time-step: +c ----------------------------------------------- + +c surface temperature + DO ig=1,ngrid + ptsoil(ig,1)=(lambda*zc(ig,1)+ptsrf(ig))/ + s (lambda*(1.-zd(ig,1))+1.) + ENDDO + +c other temperatures + DO jk=1,nsoil-1 + DO ig=1,ngrid + ptsoil(ig,jk+1)=zc(ig,jk)+zd(ig,jk)*ptsoil(ig,jk) + ENDDO + ENDDO + + ENDIF +c----------------------------------------------------------------------- +c Computation of the Cgrd and Dgrd coefficient for the next step: +c --------------------------------------------------------------- + + DO jk=1,nsoil + zdz2(jk)=dz2(jk)/ptimestep + ENDDO + + DO ig=1,ngrid + z1(ig)=zdz2(nsoil)+dz1(nsoil-1) + zc(ig,nsoil-1)=zdz2(nsoil)*ptsoil(ig,nsoil)/z1(ig) + zd(ig,nsoil-1)=dz1(nsoil-1)/z1(ig) + ENDDO + + DO jk=nsoil-1,2,-1 + DO ig=1,ngrid + z1(ig)=1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk))) + zc(ig,jk-1)= + s (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk))*z1(ig) + zd(ig,jk-1)=dz1(jk-1)*z1(ig) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c computation of the surface diffusive flux from ground and +c calorific capacity of the ground: +c --------------------------------- + + DO ig=1,ngrid + pfluxgrd(ig)=ptherm_i(ig)*dz1(1)* + s (zc(ig,1)+(zd(ig,1)-1.)*ptsoil(ig,1)) + pcapcal(ig)=ptherm_i(ig)* + s (dz2(1)+ptimestep*(1.-zd(ig,1))*dz1(1)) + z1(ig)=lambda*(1.-zd(ig,1))+1. + pcapcal(ig)=pcapcal(ig)/z1(ig) + pfluxgrd(ig)=pfluxgrd(ig) + s +pcapcal(ig)*(ptsoil(ig,1)*z1(ig)-lambda*zc(ig,1)-ptsrf(ig)) + s /ptimestep + ENDDO + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/solang.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/solang.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/solang.F (revision 4176) @@ -0,0 +1,113 @@ + subroutine solang ( kgrid,psilon,pcolon,psilat,pcolat, + & ptim1,ptim2,ptim3,pmu0,pfract ) + IMPLICIT NONE + +C +C**** *LW* - ORGANIZES THE LONGWAVE CALCULATIONS +C +C PURPOSE. +C -------- +C CALCULATES THE SOLAR ANGLE FOR ALL THE POINTS OF THE GRID +C +C** INTERFACE. +C ---------- +C SUBROUTINE SOLANG ( KGRID ) +C +C EXPLICIT ARGUMENTS : +C -------------------- +C ==== INPUTS === +C +C PSILON(KGRID) : SINUS OF THE LONGITUDE +C PCOLON(KGRID) : COSINUS OF THE LONGITUDE +C PSILAT(KGRID) : SINUS OF THE LATITUDE +C PCOLAT(KGRID) : COSINUS OF THE LATITUDE +C PTIM1 : SIN(DECLI) +C PTIM2 : COS(DECLI)*COS(TIME) +C PTIM3 : SIN(DECLI)*SIN(TIME) +C +C ==== OUTPUTS === +C +C PMU0 (KGRID) : SOLAR ANGLE +C PFRACT(KGRID) : DAY FRACTION OF THE TIME INTERVAL +C +C IMPLICIT ARGUMENTS : NONE +C -------------------- +C +C METHOD. +C ------- +C +C EXTERNALS. +C ---------- +C +C NONE +C +C REFERENCE. +C ---------- +C +C RADIATIVE PROCESSES IN METEOROLOGIE AND CLIMATOLOGIE +C PALTRIDGE AND PLATT +C +C AUTHOR. +C ------- +C FREDERIC HOURDIN +C +C MODIFICATIONS. +C -------------- +C ORIGINAL :90-01-14 +C 92-02-14 CALCULATIONS DONE THE ENTIER GRID (J.Polcher) +C----------------------------------------------------------------------- +C +C ------------------------------------------------------------------ + +C----------------------------------------------------------------------- +C +C* 0.1 ARGUMENTS +C --------- +C + INTEGER kgrid + REAL ptim1,ptim2,ptim3 + REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid) + REAL psilat(kgrid), pcolat(kgrid) +C + INTEGER jl + REAL ztim1,ztim2,ztim3 +C------------------------------------------------------------------------ +C------------------------------------------------------------------------ +C------------------------------------------------------------------------ +C +C------------------------------------------------------------------------ +C +C* 1. INITIALISATION +C -------------- +C + 100 CONTINUE +C + DO jl=1,kgrid + pmu0(jl)=0. + pfract(jl)=0. + ENDDO +C +C* 1.1 COMPUTATION OF THE SOLAR ANGLE +C ------------------------------ +C + DO jl=1,kgrid + ztim1=psilat(jl)*ptim1 + ztim2=pcolat(jl)*ptim2 + ztim3=pcolat(jl)*ptim3 + pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl) + ENDDO +C +C* 1.2 DISTINCTION BETWEEN DAY AND NIGHT +C --------------------------------- +C + DO jl=1,kgrid + IF (pmu0(jl).gt.0.) THEN + pfract(jl)=1. + ELSE + pmu0(jl)=0. + pfract(jl)=0. + ENDIF + ENDDO +C + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/solarlong.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/solarlong.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/solarlong.F (revision 4176) @@ -0,0 +1,100 @@ + SUBROUTINE solarlong(pday,psollong) + IMPLICIT NONE + +c======================================================================= +c +c Objet: +c ------ +c +c Calcul de la distance soleil-planete et de la declinaison +c en fonction du jour de l'annee. +c +c +c Methode: +c -------- +c +c Calcul complet de l'elipse +c +c Interface: +c ---------- +c +c Uncommon comprenant les parametres orbitaux. +c +c Arguments: +c ---------- +c +c Input: +c ------ +c pday jour de l'annee (le jour 0 correspondant a l'equinoxe) +c lwrite clef logique pour sorties de controle +c +c Output: +c ------- +c pdist_sol distance entre le soleil et la planete +c ( en unite astronomique pour utiliser la constante +c solaire terrestre 1370 Wm-2 ) +c pdecli declinaison ( en radians ) +c +c======================================================================= +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +#include "comorbit.h" + +c arguments: +c ---------- + + REAL pday,pdist_sol,pdecli,psollong + LOGICAL lwrite + +c Local: +c ------ + + REAL zanom,xref,zx0,zdx,zteta,zz + INTEGER iter + + +c----------------------------------------------------------------------- +c calcul de l'angle polaire et de la distance au soleil : +c ------------------------------------------------------- + +c calcul de l'zanomalie moyenne + + zz=(pday-peri_day)/year_day + zanom=2.*pi*(zz-nint(zz)) + xref=abs(zanom) + +c resolution de l'equation horaire zx0 - e * sin (zx0) = xref +c methode de Newton + + zx0=xref+e_elips*sin(xref) + DO 110 iter=1,10 + zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) + if(abs(zdx).le.(1.e-7)) goto 120 + zx0=zx0+zdx +110 continue +120 continue + zx0=zx0+zdx + if(zanom.lt.0.) zx0=-zx0 + +c zteta est la longitude solaire + + zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) + + psollong=zteta-timeperi + + IF(psollong.LT.0.) psollong=psollong+2.*pi + IF(psollong.GT.2.*pi) psollong=psollong-2.*pi +c----------------------------------------------------------------------- +c sorties eventuelles: +c --------------------- + + IF (lwrite) THEN + PRINT*,'jour de l"annee :',pday + PRINT*,'distance au soleil (en unite astronomique) :',pdist_sol + PRINT*,'declinaison (en degres) :',pdecli*180./pi + ENDIF + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/surface.h =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/surface.h (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/surface.h (revision 4176) @@ -0,0 +1,8 @@ +! REAL rnaturfi(ngridmax),phisfi(ngridmax) + REAL I_mer,I_ter,Cd_mer,Cd_ter + REAL alb_mer,alb_ter,emi_mer,emi_ter + + COMMON/surface/ +! rnaturfi,phisfi, + s I_mer,I_ter,Cd_mer,Cd_ter, + s alb_mer,alb_ter,emi_mer,emi_ter Index: /dynamico_lmdz/simple_physics/phyparam/param/sw.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/sw.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/sw.F (revision 4176) @@ -0,0 +1,227 @@ + SUBROUTINE sw(ngrid,nlayer,ldiurn, + $ coefvis,albedo, + $ plevel,ps_rad,pmu,pfract,psolarf0, + $ fsrfvis,dtsw, + $ lwrite) + IMPLICIT NONE +c======================================================================= +c +c Rayonnement solaire en atmosphere non diffusante avec un +c coefficient d'absoprption gris. +c +c======================================================================= +c +c declarations: +c ------------- +c +#include "comcstfi.h" +c +c arguments: +c ---------- +c + INTEGER ngrid,nlayer + REAL albedo(ngrid),coefvis + REAL pmu(ngrid),pfract(ngrid) + REAL plevel(ngrid,nlayer+1),ps_rad + REAL psolarf0 + REAL fsrfvis(ngrid),dtsw(ngrid,nlayer) + LOGICAL lwrite,ldiurn +c +c variables locales: +c ------------------ +c + + REAL zalb(ngrid),zmu(ngrid),zfract(ngrid) + REAL zplev(ngrid,nlayer+1) + REAL zflux(ngrid),zdtsw(ngrid,nlayer) + + INTEGER ig,l,nlevel,index(ngrid),ncount,igout + REAL ztrdir(ngrid,nlayer+1),ztrref(ngrid,nlayer+1) + REAL zfsrfref(ngrid) + REAL z1(ngrid) + REAL zu(ngrid,nlayer+1) + REAL tau0 + + LOGICAL firstcall + SAVE firstcall + DATA firstcall/.true./ +!$OMP THREADPRIVATE(firstcall) + +c----------------------------------------------------------------------- +c 1. initialisations: +c ------------------- + + + nlevel=nlayer+1 + +c----------------------------------------------------------------------- +c Definitions des tableaux locaux pour les points ensoleilles: +c ------------------------------------------------------------ + + IF (ldiurn) THEN + ncount=0 + DO ig=1,ngrid + index(ig)=0 + ENDDO + DO ig=1,ngrid + IF(pfract(ig).GT.1.e-6) THEN + ncount=ncount+1 + index(ncount)=ig + ENDIF + ENDDO + CALL monGATHER(ncount,zfract,pfract,index) + CALL monGATHER(ncount,zmu,pmu,index) + CALL monGATHER(ncount,zalb,albedo,index) + DO l=1,nlevel + CALL monGATHER(ncount,zplev(1,l),plevel(1,l),index) + ENDDO + ELSE + ncount=ngrid + zfract(:)=pfract(:) + zmu(:)=pmu(:) + zalb(:)=albedo(:) + zplev(:,:)=plevel(:,:) + ENDIF + +c----------------------------------------------------------------------- +c calcul des profondeurs optiques integres depuis p=0: +c ---------------------------------------------------- + + tau0=-.5*log(coefvis) + +c calcul de la partie homogene de l'opacite + tau0=tau0/ps_rad + DO l=1,nlayer+1 + DO ig=1,ncount + zu(ig,l)=tau0*zplev(ig,l) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 2. calcul de la transmission depuis le sommet de l'atmosphere: +c ----------------------------------------------------------- + + DO l=1,nlevel + DO ig=1,ncount + ztrdir(ig,l)=exp(-zu(ig,l)/zmu(ig)) + ENDDO + ENDDO + + IF (lwrite) THEN + igout=ncount/2+1 + PRINT* + PRINT*,'Diagnostique des transmission dans le spectre solaire' + PRINT*,'zfract, zmu, zalb' + PRINT*,zfract(igout), zmu(igout), zalb(igout) + PRINT*,'Pression, quantite d abs, transmission' + DO l=1,nlayer+1 + PRINT*,zplev(igout,l),zu(igout,l),ztrdir(igout,l) + ENDDO + ENDIF + +c----------------------------------------------------------------------- +c 3. taux de chauffage, ray. solaire direct: +c ------------------------------------------ + + DO l=1,nlayer + DO ig=1,ncount + zdtsw(ig,l)=g*psolarf0*zfract(ig)*zmu(ig)* + $ (ztrdir(ig,l+1)-ztrdir(ig,l))/ + $ (cpp*(zplev(ig,l)-zplev(ig,l+1))) + ENDDO + ENDDO + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires:' + PRINT*,' 1 taux de chauffage lie au ray. solaire direct' + DO l=1,nlayer + PRINT*,zdtsw(igout,l) + ENDDO + ENDIF + + +c----------------------------------------------------------------------- +c 4. calcul du flux solaire arrivant sur le sol: +c ---------------------------------------------- + + DO ig=1,ncount + z1(ig)=zfract(ig)*zmu(ig)*psolarf0*ztrdir(ig,1) + zflux(ig)=(1.-zalb(ig))*z1(ig) + zfsrfref(ig)= zalb(ig)*z1(ig) + ENDDO + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires:' + PRINT*,' 2 flux solaire net incident sur le sol' + PRINT*,zflux(igout) + ENDIF + + +c----------------------------------------------------------------------- +c 5.calcul des traansmissions depuis le sol, cas diffus: +c ------------------------------------------------------ + + DO l=1,nlevel + DO ig=1,ncount + ztrref(ig,l)=exp(-(zu(ig,1)-zu(ig,l))*1.66) + ENDDO + ENDDO + + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires' + PRINT*,' 3 transmission avec les sol' + PRINT*,'niveau transmission' + DO l=1,nlevel + PRINT*,l,ztrref(igout,l) + ENDDO + ENDIF + +c----------------------------------------------------------------------- +c 6.ajout a l'echauffement de la contribution du ray. sol. reflechit: +c ------------------------------------------------------------------- + + DO l=1,nlayer + DO ig=1,ncount + zdtsw(ig,l)=zdtsw(ig,l)+ + $ g*zfsrfref(ig)*(ztrref(ig,l+1)-ztrref(ig,l))/ + $ (cpp*(zplev(ig,l+1)-zplev(ig,l))) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 10. sorites eventuelles: +c ------------------------ + + IF (lwrite) THEN + PRINT* + PRINT*,'Diagnostique des taux de chauffage solaires:' + PRINT*,' 3 taux de chauffage total' + DO l=1,nlayer + PRINT*,zdtsw(igout,l) + ENDDO + ENDIF + + IF (ldiurn) THEN + CALL zerophys(ngrid,fsrfvis) + CALL monscatter(ncount,fsrfvis,index,zflux) + CALL zerophys(ngrid*nlayer,dtsw) + DO l=1,nlayer + CALL monscatter(ncount,dtsw(1,l),index,zdtsw(1,l)) + ENDDO + ELSE + print*,'NOT DIURNE' + fsrfvis(:)=zflux(:) + dtsw(:,:)=zdtsw(:,:) + ENDIF +c call dump2d(iim,jjm-1,zflux(2),'ZFLUX ') +c call dump2d(iim,jjm-1,fsrfvis(2),'FSRVIS ') +c call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') +c call dump2d(iim,jjm-1,pmu(2),'pmu ') +c call dump2d(iim,jjm-1,pfract(2),'pfract ') +c call dump2d(iim,jjm-1,albedo(2),'albedo ') +c call dump2d(iim,jjm-1,ztrdir(2,1),'ztrdir ') + + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/vdif.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/vdif.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/vdif.F (revision 4176) @@ -0,0 +1,305 @@ + subroutiNE vdif(ngrid,nlay,ptime, + $ ptimestep,pcapcal,pz0, + $ pplay,pplev,pzlay,pzlev, + $ pu,pv,ph,ptsrf,pemis, + $ pdufi,pdvfi,pdhfi,pfluxsrf, + $ pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, + $ lwrite) + IMPLICIT NONE + +c======================================================================= +c +c Diffusion verticale +c Shema implicite +c On commence par rajouter au variables x la tendance physique +c et on resoult en fait: +c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) +c +c arguments: +c ---------- +c +c entree: +c ------- +c +c +c======================================================================= + +c----------------------------------------------------------------------- +c declarations: +c ------------- + +#include "comcstfi.h" +c des include de la geometrie dynamique pour sorties graphiques +! #include "paramet.h" +! #include "comvert.h" +! #include "comgeom.h" +#include "description.h" +c +c arguments: +c ---------- + + INTEGER ngrid,nlay + REAL ptime,ptimestep + REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) + REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1) + REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay) + REAL ptsrf(ngrid),pemis(ngrid) + REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay) + REAL pfluxsrf(ngrid) + REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay) + REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid) + REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1) + LOGICAL lwrite +c +c local: +c ------ + + INTEGER ilev,ig,ilay,nlev + INTEGER unit,ierr,it1,it2 + INTEGER cluvdb,putdat,putvdim,setname,setvdim + REAL z4st,zdplanck(ngrid),zu2 + REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1) + REAL zcdv(ngrid),zcdh(ngrid) + REAL zu(ngrid,nlay),zv(ngrid,nlay) + REAL zh(ngrid,nlay) + REAL ztsrf2(ngrid) + REAL z1(ngrid),z2(ngrid) + REAL za(ngrid,nlay),zb(ngrid,nlay) + REAL zb0(ngrid,nlay) + REAL zc(ngrid,nlay),zd(ngrid,nlay) + REAL zcst1 + + EXTERNAL coefdifv + REAL, SAVE :: karman=0.4 + INTEGER, SAVE :: icount=0 +!$OMP THREADPRIVATE(karman,icount) + +c +c----------------------------------------------------------------------- +c initialisations: +c ---------------- + + nlev=nlay+1 + +c computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp: +c with rho=p/RT=p/ (R Theta) (p/ps)**kappa +c --------------------------------- + + DO ilay=1,nlay + DO ig=1,ngrid + za(ig,ilay)= + s (pplev(ig,ilay)-pplev(ig,ilay+1))/g + ENDDO + ENDDO + + zcst1=4.*g*ptimestep/(r*r) + DO ilev=2,nlev-1 + DO ig=1,ngrid + zb0(ig,ilev)=pplev(ig,ilev)* + s (pplev(ig,1)/pplev(ig,ilev))**rcp / + s (ph(ig,ilev-1)+ph(ig,ilev)) + zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/ + s (pplay(ig,ilev-1)-pplay(ig,ilev)) + ENDDO + ENDDO + DO ig=1,ngrid + zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) + ENDDO + IF(lwrite) THEN + ig=ngrid/2+1 + PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz' + DO ilay=1,nlay + WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay), + s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay) + ENDDO + PRINT*,'Pression (mbar) ,altitude (km),zb' + DO ilev=1,nlay + WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev), + s zb0(ig,ilev) + ENDDO + ENDIF + +c----------------------------------------------------------------------- +c 2. ajout des tendances physiques: +c ------------------------------ + + DO ilev=1,nlay + DO ig=1,ngrid + zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep + zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep + zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c 3. calcul de cd : +c ---------------- +c + CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh) +c CALL my_25(ngrid,nlay,ptimestep,g,pzlev,pzlay,pu,pv,ph,zcdv, +c a pq2,pq2l,zkv,zkh) + CALL vdif_k(ngrid,nlay, + s ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zcdv,zkv,zkh,pq2,pq2l) + + DO ig=1,ngrid + zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) + zcdv(ig)=zcdv(ig)*sqrt(zu2) + zcdh(ig)=zcdh(ig)*sqrt(zu2) + ENDDO + + IF(lwrite) THEN + PRINT* + PRINT*,'Diagnostique diffusion verticale' + PRINT*,'coefficients Cd pour v et h' + PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1) + PRINT*,'coefficients K pour v et h' + DO ilev=1,nlay + PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev) + ENDDO + ENDIF + +c----------------------------------------------------------------------- +c integration verticale pour u: +c ----------------------------- +c + CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2)) + CALL multipl(ngrid,zcdv,zb0,zb) + DO ig=1,ngrid + z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) + zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig) + zd(ig,nlay)=zb(ig,nlay)*z1(ig) + ENDDO + + DO ilay=nlay-1,1,-1 + DO ig=1,ngrid + z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ + $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) + zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ + $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) + zd(ig,ilay)=zb(ig,ilay)*z1(ig) + ENDDO + ENDDO + + DO ig=1,ngrid + zu(ig,1)=zc(ig,1) + ENDDO + DO ilay=2,nlay + DO ig=1,ngrid + zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c integration verticale pour v: +c ----------------------------- +c + DO ig=1,ngrid + z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) + zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig) + zd(ig,nlay)=zb(ig,nlay)*z1(ig) + ENDDO + + DO ilay=nlay-1,1,-1 + DO ig=1,ngrid + z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ + $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) + zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ + $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) + zd(ig,ilay)=zb(ig,ilay)*z1(ig) + ENDDO + ENDDO + + DO ig=1,ngrid + zv(ig,1)=zc(ig,1) + ENDDO + DO ilay=2,nlay + DO ig=1,ngrid + zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c integration verticale pour h: +c ----------------------------- +c + CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) + CALL multipl(ngrid,zcdh,zb0,zb) + DO ig=1,ngrid + z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) + zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig) + zd(ig,nlay)=zb(ig,nlay)*z1(ig) + ENDDO + + DO ilay=nlay-1,1,-1 + DO ig=1,ngrid + z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ + $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) + zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ + $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) + zd(ig,ilay)=zb(ig,ilay)*z1(ig) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c rajout eventuel de planck dans le shema implicite: +c -------------------------------------------------- + + z4st=4.*5.67e-8*ptimestep +c z4st=0. + DO ig=1,ngrid + zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) + ENDDO + +c----------------------------------------------------------------------- +c calcul le l'evolution de la temperature du sol: +c ----------------------------------------------- + + DO ig=1,ngrid + z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) + s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep + z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) + ztsrf2(ig)=z1(ig)/z2(ig) + zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) + pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep + ENDDO + +c----------------------------------------------------------------------- +c integration verticale finale: +c ----------------------------- + + DO ilay=2,nlay + DO ig=1,ngrid + zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1) + ENDDO + ENDDO + +c----------------------------------------------------------------------- +c calcul final des tendances de la diffusion verticale: +c ----------------------------------------------------- + + DO ilev = 1, nlay + DO ig=1,ngrid + pdudif(ig,ilev)=( zu(ig,ilev)- + $ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep + pdvdif(ig,ilev)=( zv(ig,ilev)- + $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep + pdhdif(ig,ilev)=( zh(ig,ilev)- + $ (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep + ENDDO + ENDDO + + IF(lwrite) THEN + PRINT* + PRINT*,'Diagnostique de la diffusion verticale' + PRINT*,'h avant et apres diffusion verticale' + PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1) + DO 3110 ilev=1,nlay + PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev) +3110 CONTINUE + ENDIF + + +c----------------------------------------------------------------------- + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/vdif_cd.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/vdif_cd.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/vdif_cd.F (revision 4176) @@ -0,0 +1,97 @@ + SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh) + IMPLICIT NONE +c======================================================================= +c +c Subject: computation of the surface drag coefficient using the +c ------- approch developed by Loui for ECMWF. +c +c Author: Frederic Hourdin 15 /10 /93 +c ------- +c +c Arguments: +c ---------- +c +c inputs: +c ------ +c ngrid size of the horizontal grid +c pg gravity (m s -2) +c pz(ngrid) height of the first atmospheric layer +c pu(ngrid) u component of the wind in that layer +c pv(ngrid) v component of the wind in that layer +c pts(ngrid) surfacte temperature +c ph(ngrid) potential temperature T*(p/ps)^kappa +c +c outputs: +c -------- +c pcdv(ngrid) Cd for the wind +c pcdh(ngrid) Cd for potential temperature +c +c======================================================================= +c +c----------------------------------------------------------------------- +c Declarations: +c ------------- + +c Arguments: +c ---------- + + INTEGER ngrid,nlay + REAL pz0(ngrid) + REAL pg,pz(ngrid) + REAL pu(ngrid),pv(ngrid) + REAL pts(ngrid),ph(ngrid) + REAL pcdv(ngrid),pcdh(ngrid) + +c Local: +c ------ + + INTEGER ig + + REAL zu2,z1,zri,zcd0,zz + + REAL karman,b,c,d,c2b,c3bc,c3b,z0,umin2 + LOGICAL firstcal + DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/ + DATA firstcal/.true./ + SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2 + +c----------------------------------------------------------------------- +c couche de surface: +c ------------------ + +c DO ig=1,ngrid +c zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2 +c pcdv(ig)=pz0(ig)*(1.+sqrt(zu2)) +c pcdh(ig)=pcdv(ig) +c ENDDO +c RETURN + + IF (firstcal) THEN + c2b=2.*b + c3bc=3.*b*c + c3b=3.*b + firstcal=.false. + ENDIF + +c!!!! WARNING, verifier la formule originale de Louis! + DO ig=1,ngrid + zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2 + zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2) + z1=1.+pz(ig)/pz0(ig) + zcd0=karman/log(z1) + zcd0=zcd0*zcd0*sqrt(zu2) + IF(zri.LT.0.) THEN + z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri)) + pcdv(ig)=zcd0*(1.-2.*z1) + pcdh(ig)=zcd0*(1.-3.*z1) + ELSE + zz=sqrt(1.+d*zri) + pcdv(ig)=zcd0/(1.+c2b*zri/zz) + pcdh(ig)=zcd0/(1.+c3b*zri*zz) + ENDIF + ENDDO + +c----------------------------------------------------------------------- + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/vdif_k.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/vdif_k.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/vdif_k.F (revision 4176) @@ -0,0 +1,59 @@ + SUBROUTINE vdif_k(ngrid,nlay, + s ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pcdv,pkv,pkh) + + IMPLICIT NONE +#include "planete.h" + + INTEGER ngrid,nlay + + REAL ptimestep + REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1) + REAL pz0(ngrid) + REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay) + REAL pg,pcdv(ngrid) + REAL pkv(ngrid,nlay+1),pkh(ngrid,nlay+1) + + INTEGER ig,il + REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix + + REAL karman + SAVE karman + DATA lmixmin,emin_turb,karman/100.,1.e-8,.4/ + + + print*,'LMIXMIN',lmixmin + DO ig=1,ngrid + pkv(ig,1)=0. + pkh(ig,1)=0. + pkv(ig,nlay+1)=0. + pkh(ig,nlay+1)=0. + ENDDO +c s ' zdu,zdv,zdz,zdovdz2,ph(ig,il)+ph(ig,il-1)' + DO il=2,nlay + DO ig=1,ngrid + z1=pzlev(ig,il)+pz0(ig) + lmix=karman*z1/(1.+karman*z1/lmixmin) +c lmix=lmixmin +c WARNING test lmix=lmixmin + zdu=pu(ig,il)-pu(ig,il-1) + zdv=pv(ig,il)-pv(ig,il-1) + zdz=pzlay(ig,il)-pzlay(ig,il-1) + zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz) + IF(zdvodz2.LT.1.e-5) THEN + pkv(ig,il)=lmix*sqrt(emin_turb) + ELSE + zri=2.*pg*(ph(ig,il)-ph(ig,il-1)) + s / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2 ) + pkv(ig,il)= + s lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb)) + ENDIF + pkh(ig,il)=pkv(ig,il) +c IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il), +c s zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1), +c s lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1), +c s ph(ig,il),ph(ig,il-1) + ENDDO + ENDDO + + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/param/zenang.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/zenang.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/zenang.F (revision 4176) @@ -0,0 +1,219 @@ + SUBROUTINE zenang(klon,longi,gmtime,pdtrad,lat,long, + s pmu0,frac) + IMPLICIT none +c============================================================= +c Auteur : O. Boucher (LMD/CNRS) +c d'apres les routines zenith et angle de Z.X. Li +c Objet : calculer les valeurs moyennes du cos de l'angle zenithal +c et l'ensoleillement moyen entre gmtime1 et gmtime2 +c connaissant la declinaison, la latitude et la longitude. +c Rque : Different de la routine angle en ce sens que zenang +c fournit des moyennes de pmu0 et non des valeurs +c instantanees, du coup frac prend toutes les valeurs +c entre 0 et 1. +c Date : premiere version le 13 decembre 1994 +c revu pour GCM le 30 septembre 1996 +c=============================================================== +c longi----INPUT : la longitude vraie de la terre dans son plan +c solaire a partir de l'equinoxe de printemps (degre) +c gmtime---INPUT : temps universel en fraction de jour +c pdtrad---INPUT : pas de temps du rayonnement (secondes) +c lat------INPUT : latitude en degres +c long-----INPUT : longitude en degres +c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad +c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad +c================================================================ +#include "comorbit.h" + integer klon +c================================================================ + real longi, gmtime, pdtrad + real lat(klon), long(klon), pmu0(klon), frac(klon) +c================================================================ + integer i + real gmtime1, gmtime2 + real pi_local, deux_pi_local, incl + real omega1, omega2, omega +c omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. +c omega : heure en radian du coucher de soleil +c -omega est donc l'heure en radian de lever du soleil + real omegadeb, omegafin + real zfrac1, zfrac2, z1_mu, z2_mu + real lat_sun ! declinaison en radian + real lon_sun ! longitude solaire en radian + real latr ! latitude du pt de grille en radian +c================================================================ +c + pi_local = 4.0 * ATAN(1.0) + deux_pi_local = 2.0 * pi_local +c incl=R_incl * pi_local / 180. + print*,'Obliquite =' ,obliquit + incl=obliquit * pi_local / 180. +c +c lon_sun = longi * pi_local / 180.0 + lon_sun = longi + lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) +c + gmtime1=gmtime*86400. + gmtime2=gmtime*86400.+pdtrad +c + DO i = 1, klon +c +c latr = lat(i) * pi_local / 180. + latr = lat(i) +c +c--pose probleme quand lat=+/-90 degres +c +c omega = -TAN(latr)*TAN(lat_sun) +c omega = ACOS(omega) +c IF (latr.GE.(pi_local/2.+lat_sun) +c . .OR. latr.LE.(-pi_local/2.+lat_sun)) THEN +c omega = 0.0 ! nuit polaire +c ENDIF +c IF (latr.GE.(pi_local/2.-lat_sun) +c . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN +c omega = pi_local ! journee polaire +c ENDIF +c +c--remplace par cela (le cas par defaut est different) +c + omega=0.0 !--nuit polaire + IF (latr.GE.(pi_local/2.-lat_sun) + . .OR. latr.LE.(-pi_local/2.-lat_sun)) THEN + omega = pi_local ! journee polaire + ENDIF + IF (latr.LT.(pi_local/2.+lat_sun).AND. + . latr.GT.(-pi_local/2.+lat_sun).AND. + . latr.LT.(pi_local/2.-lat_sun).AND. + . latr.GT.(-pi_local/2.-lat_sun)) THEN + omega = -TAN(latr)*TAN(lat_sun) + omega = ACOS(omega) + ENDIF +c + omega1 = gmtime1 + long(i)*86400.0/360.0 + omega1 = omega1 / 86400.0*deux_pi_local + omega1 = MOD (omega1+deux_pi_local, deux_pi_local) + omega1 = omega1 - pi_local +c + omega2 = gmtime2 + long(i)*86400.0/360.0 + omega2 = omega2 / 86400.0*deux_pi_local + omega2 = MOD (omega2+deux_pi_local, deux_pi_local) + omega2 = omega2 - pi_local +c + IF (omega1.LE.omega2) THEN !--on est dans la meme journee locale +c + IF (omega2.LE.-omega .OR. omega1.GE.omega + . .OR. omega.LT.1e-5) THEN !--nuit + frac(i)=0.0 + pmu0(i)=0.0 + ELSE !--jour+nuit/jour + omegadeb=MAX(-omega,omega1) + omegafin=MIN(omega,omega2) + frac(i)=(omegafin-omegadeb)/(omega2-omega1) + pmu0(i)=SIN(latr)*SIN(lat_sun) + + . COS(latr)*COS(lat_sun)* + . (SIN(omegafin)-SIN(omegadeb))/ + . (omegafin-omegadeb) + ENDIF +c + ELSE !---omega1 GT omega2 -- a cheval sur deux journees +c +c-------------------entre omega1 et pi + IF (omega1.GE.omega) THEN !--nuit + zfrac1=0.0 + z1_mu =0.0 + ELSE !--jour+nuit + omegadeb=MAX(-omega,omega1) + omegafin=omega + zfrac1=omegafin-omegadeb + z1_mu =SIN(latr)*SIN(lat_sun) + + . COS(latr)*COS(lat_sun)* + . (SIN(omegafin)-SIN(omegadeb))/ + . (omegafin-omegadeb) + ENDIF +c---------------------entre -pi et omega2 + IF (omega2.LE.-omega) THEN !--nuit + zfrac2=0.0 + z2_mu =0.0 + ELSE !--jour+nuit + omegadeb=-omega + omegafin=MIN(omega,omega2) + zfrac2=omegafin-omegadeb + z2_mu =SIN(latr)*SIN(lat_sun) + + . COS(latr)*COS(lat_sun)* + . (SIN(omegafin)-SIN(omegadeb))/ + . (omegafin-omegadeb) +c + ENDIF +c-----------------------moyenne + frac(i)=(zfrac1+zfrac2)/(omega2+deux_pi_local-omega1) + pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10) +c + ENDIF !---comparaison omega1 et omega2 +c + ENDDO +c + END +c=================================================================== +#ifdef PASLA + SUBROUTINE zenith (klon,longi, gmtime, lat, long, + s pmu0, fract) + IMPLICIT none +c +c Auteur(s): Z.X. Li (LMD/ENS) +c +c Objet: calculer le cosinus de l'angle zenithal du soleil en +c connaissant la declinaison du soleil, la latitude et la +c longitude du point sur la terre, et le temps universel +c +c Arguments d'entree: +c longi : declinaison du soleil (en degres) +c gmtime : temps universel en second qui varie entre 0 et 86400 +c lat : latitude en degres +c long : longitude en degres +c Arguments de sortie: +c pmu0 : cosinus de l'angle zenithal +c +c==================================================================== +#include "YOMCST.h" +c==================================================================== + INTETER klon + REAL longi, gmtime + REAL lat(klon), long(klon), pmu0(klon), fract(klon) +c===================================================================== + INTEGER n + REAL zpi, zpir, omega, zgmtime + REAL incl, lat_sun, lon_sun +c---------------------------------------------------------------------- + zpi = 4.0*ATAN(1.0) + zpir = zpi / 180.0 + zgmtime=gmtime*86400. +c + incl=R_incl * zpir +c + lon_sun = longi * zpir + lat_sun = ASIN (SIN(lon_sun)*SIN(incl) ) +c +c--initialisation a la nuit +c + DO n =1, klon + pmu0(n)=0. + fract(n)=0.0 + ENDDO +c +c 1 degre en longitude = 240 secondes en temps +c + DO n = 1, klon + omega = zgmtime + long(n)*86400.0/360.0 + omega = omega / 86400.0 * 2.0 * zpi + omega = MOD(omega + 2.0 * zpi, 2.0 * zpi) + omega = omega - zpi + pmu0(n) = sin(lat(n)*zpir) * sin(lat_sun) + . + cos(lat(n)*zpir) * cos(lat_sun) + . * cos(omega) + pmu0(n) = MAX (pmu0(n), 0.0) + IF (pmu0(n).GT.1.E-6) fract(n)=1.0 + ENDDO +c + RETURN + END +#endif Index: /dynamico_lmdz/simple_physics/phyparam/param/zerofi.F =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/param/zerofi.F (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/param/zerofi.F (revision 4176) @@ -0,0 +1,11 @@ + SUBROUTINE zerophys(n,x) + IMPLICIT NONE + INTEGER n + REAL x(n) + INTEGER i + + DO i=1,n + x(i)=0. + ENDDO + RETURN + END Index: /dynamico_lmdz/simple_physics/phyparam/physiq_mod.F90 =================================================================== --- /dynamico_lmdz/simple_physics/phyparam/physiq_mod.F90 (revision 4176) +++ /dynamico_lmdz/simple_physics/phyparam/physiq_mod.F90 (revision 4176) @@ -0,0 +1,224 @@ +! $Id: physiq.F 1565 2011-08-31 12:53:29Z jghattas $ +MODULE physiq_mod + +IMPLICIT NONE + +CONTAINS + + SUBROUTINE physiq (nlon,nlev, & + & debut,lafin,pdtphys, & + & paprs,pplay,pphi,pphis,presnivs, & + & u,v,t,qx, & + & flxmass_w, & + & d_u, d_v, d_t, d_qx, d_ps) + + USE dimphy, only : klon,klev + USE infotrac_phy, only : nqtot + USE geometry_mod, only : latitude + USE comcstphy, only : rg + USE iophy, only : histbeg_phy,histwrite_phy + USE ioipsl, only : getin,histvert,histdef,histend,ymds2ju + USE mod_phys_lmdz_para, only : jj_nb + USE phys_state_var_mod, only : phys_state_var_init + USE mod_grid_phy_lmdz, ONLY: nbp_lon,nbp_lat + +#ifdef CPP_XIOS + USE xios, ONLY: xios_update_calendar + USE wxios, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef + USE iophy, ONLY: histwrite_phy +#endif + + IMPLICIT none +! +! Routine argument: +! + integer,intent(in) :: nlon ! number of atmospheric colums + integer,intent(in) :: nlev ! number of vertical levels (should be =klev) + logical,intent(in) :: debut ! signals first call to physics + logical,intent(in) :: lafin ! signals last call to physics + real,intent(in) :: pdtphys ! physics time step (s) + real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) + real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) + real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer + real,intent(in) :: pphis(klon) ! surface geopotential + real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers + real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s) + real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s) + real,intent(in) :: t(klon,klev) ! temperature (K) + real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air) + real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux + real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s) + real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s) + real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s) + real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers + real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure + +integer,save :: itau=0 ! counter to count number of calls to physics +!$OMP THREADPRIVATE(itau) +real :: temp_newton(klon,klev) +integer :: k +logical, save :: first=.true. +!$OMP THREADPRIVATE(first) + +! For I/Os +integer :: itau0 +real :: zjulian +real :: dtime +integer :: nhori ! horizontal coordinate ID +integer,save :: nid_hist ! output file ID +!$OMP THREADPRIVATE(nid_hist) +integer :: zvertid ! vertical coordinate ID +integer,save :: iwrite_phys ! output every iwrite_phys physics step +!$OMP THREADPRIVATE(iwrite_phys) +integer,save :: iwrite_phys_omp ! intermediate variable to read iwrite_phys + ! (must be shared by all threads) +real :: t_ops ! frequency of the IOIPSL operations (eg average over...) +real :: t_wrt ! frequency of the IOIPSL outputs +REAL,SAVE :: rjourvrai=0. +REAL,SAVE :: gmtime=0. +!$OMP THREADPRIVATE(rjourvrai,gmtime) + + +! initializations +if (debut) then ! Things to do only for the first call to physics +! load initial conditions for physics (including the grid) + call phys_state_var_init() ! some initializations, required before calling phyetat0 + call phyetat0("startphy.nc") + +! Initialize outputs: + itau0=0 +!$OMP MASTER + iwrite_phys_omp=1 !default: output every physics timestep + ! NB: getin() is not threadsafe; only one thread should call it. + call getin("iwrite_phys",iwrite_phys_omp) +!$OMP END MASTER +!$OMP BARRIER + iwrite_phys=iwrite_phys_omp + t_ops=pdtphys*iwrite_phys ! frequency of the IOIPSL operation + t_wrt=pdtphys*iwrite_phys ! frequency of the outputs in the file + ! compute zjulian for annee0=1979 and month=1 dayref=1 and hour=0.0 + !CALL ymds2ju(annee0, month, dayref, hour, zjulian) + call ymds2ju(1979, 1, 1, 0.0, zjulian) + dtime=pdtphys +#ifndef CPP_IOIPSL_NO_OUTPUT + ! Initialize IOIPSL output file + call histbeg_phy("histins.nc",itau0,zjulian,dtime,nhori,nid_hist) +#endif + +!$OMP MASTER + +#ifndef CPP_IOIPSL_NO_OUTPUT +! IOIPSL + ! define vertical coordinate + call histvert(nid_hist,"presnivs","Vertical levels","Pa",klev, & + presnivs,zvertid,'down') + ! define variables which will be written in "histins.nc" file + call histdef(nid_hist,'temperature','Atmospheric temperature','K', & + nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & + 'inst(X)',t_ops,t_wrt) + call histdef(nid_hist,'u','Eastward Zonal Wind','m/s', & + nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & + 'inst(X)',t_ops,t_wrt) + call histdef(nid_hist,'v','Northward Meridional Wind','m/s', & + nbp_lon,jj_nb,nhori,klev,1,klev,zvertid,32, & + 'inst(X)',t_ops,t_wrt) + call histdef(nid_hist,'ps','Surface Pressure','Pa', & + nbp_lon,jj_nb,nhori,1,1,1,zvertid,32, & + 'inst(X)',t_ops,t_wrt) + ! end definition sequence + call histend(nid_hist) +#endif + +#ifdef CPP_XIOS +!XIOS + ! Declare available vertical axes to be used in output files: + CALL wxios_add_vaxis("presnivs", klev, presnivs) + + ! Declare calendar and time step + CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0) + + !Finalize the context: + CALL wxios_closedef() +#endif +!$OMP END MASTER +!$OMP BARRIER +endif ! of if (debut) + +! increment local time counter itau +itau=itau+1 +gmtime=gmtime+pdtphys/86400. +IF (gmtime>=1.) THEN + gmtime=gmtime-1. + rjourvrai=rjourvrai+1. +ENDIF + +IF ( 1 == 0 ) THEN +! set all tendencies to zero +d_u(1:klon,1:klev)=0. +d_v(1:klon,1:klev)=0. +d_t(1:klon,1:klev)=0. +d_qx(1:klon,1:klev,1:nqtot)=0. +d_ps(1:klon)=0. + +! compute tendencies to return to the dynamics: +! "friction" on the first layer +d_u(1:klon,1)=-u(1:klon,1)/86400. +d_v(1:klon,1)=-v(1:klon,1)/86400. +! newtonian relaxation towards temp_newton() +do k=1,klev + temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 + d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5 +enddo + +ELSE + + CALL phyparam(klon,klev,nqtot, & + debut,lafin, & + rjourvrai,gmtime,pdtphys, & + paprs,pplay,pphi,pphis,presnivs, & + u,v,t,qx, & + flxmass_w, & + d_u,d_v,d_t,d_qx,d_ps) + + +ENDIF + + +print*,'PHYDEV: itau=',itau + +! write some outputs: +! IOIPSL +#ifndef CPP_IOIPSL_NO_OUTPUT +if (modulo(itau,iwrite_phys)==0) then + call histwrite_phy(nid_hist,.false.,"temperature",itau,t) + call histwrite_phy(nid_hist,.false.,"u",itau,u) + call histwrite_phy(nid_hist,.false.,"v",itau,v) + call histwrite_phy(nid_hist,.false.,"ps",itau,paprs(:,1)) +endif +#endif + +!XIOS +#ifdef CPP_XIOS +!$OMP MASTER + !Increment XIOS time + CALL xios_update_calendar(itau) +!$OMP END MASTER +!$OMP BARRIER + + !Send fields to XIOS: (NB these fields must also be defined as + ! in iodef.xml to be correctly used + CALL histwrite_phy("temperature",t) + CALL histwrite_phy("temp_newton",temp_newton) + CALL histwrite_phy("u",u) + CALL histwrite_phy("v",v) + CALL histwrite_phy("ps",paprs(:,1)) +#endif + +! if lastcall, then it is time to write "restartphy.nc" file +if (lafin) then + call phyredem("restartphy.nc") +endif + +end subroutine physiq + +END MODULE physiq_mod