Index: LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F (revision 1748) +++ LMDZ5/trunk/libf/dyn3dmem/calfis_loc.F (revision 1749) @@ -34,8 +34,8 @@ USE dimphy USE mod_phys_lmdz_para, mpi_root_xx=>mpi_root - USE parallel, ONLY : omp_chunk, using_mpi,jjb_u,jje_u,jjb_v,jje_v USE mod_interface_dyn_phys USE IOPHY #endif + USE parallel, ONLY : omp_chunk, using_mpi,jjb_u,jje_u,jjb_v,jje_v USE Write_Field Use Write_field_p @@ -116,6 +116,6 @@ c ----------- LOGICAL lafin - REAL heure - +! REAL heure + REAL, intent(in):: jD_cur, jH_cur REAL pvcov(iip1,jjb_v:jje_v,llm) REAL pucov(iip1,jjb_u:jje_u,llm) @@ -130,4 +130,5 @@ REAL pdteta(iip1,jjb_u:jje_u,llm) REAL pdq(iip1,jjb_u:jje_u,llm,nqtot) + REAL flxw(iip1,jjb_u:jje_u,llm) ! Flux de masse verticale sur la grille dynamique c REAL pps(iip1,jjb_u:jje_u) @@ -226,5 +227,4 @@ REAL PVteta(klon,ntetaSTD) - REAL flxw(iip1,jjb_u:jje_u,llm) ! Flux de masse verticale sur la grille dynamique REAL SSUM @@ -234,5 +234,4 @@ SAVE firstcal,debut c$OMP THREADPRIVATE(firstcal,debut) - REAL, intent(in):: jD_cur, jH_cur REAL,SAVE,dimension(1:iim,1:llm):: du_send,du_recv,dv_send,dv_recv Index: LMDZ5/trunk/libf/dyn3dmem/divgrad_p.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/divgrad_p.F (revision 1748) +++ (revision ) @@ -1,91 +1,0 @@ - SUBROUTINE divgrad_p (klevel,h, lh, divgra_out ) - USE parallel - USE times - IMPLICIT NONE -c -c======================================================================= -c -c Auteur : P. Le Van -c ---------- -c -c lh -c calcul de (div( grad )) de h ..... -c h et lh sont des arguments d'entree pour le s-prog -c divgra est un argument de sortie pour le s-prog -c -c======================================================================= -c -c declarations: -c ------------- -c -#include "dimensions.h" -#include "paramet.h" -#include "comgeom.h" -#include "comdissipn.h" -#include "logic.h" -c - INTEGER klevel - REAL h( ip1jmp1,klevel ), divgra_out( ip1jmp1,klevel ) - REAL,SAVE :: divgra( ip1jmp1,llm ) - -c - REAL ghy(ip1jm,llm), ghx(ip1jmp1,llm) - - INTEGER l,ij,iter,lh -c - INTEGER ijb,ije,jjb,jje -c -c -c CALL SCOPY ( ip1jmp1*klevel,h,1,divgra,1 ) - - ijb=ij_begin - ije=ij_end -c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) - DO l = 1, klevel - divgra(ijb:ije,l)=h(ijb:ije,l) - ENDDO -c$OMP END DO NOWAIT -c - -c - DO 10 iter = 1,lh - - jjb=jj_begin - jje=jj_end - CALL filtreg_p ( divgra,jjb,jje,jjp1,klevel,2,1,.true.,1 ) - -c call exchange_Hallo(divgra,ip1jmp1,llm,0,1) -c$OMP BARRIER -c$OMP MASTER - call suspend_timer(timer_dissip) - call exchange_Hallo(divgra,ip1jmp1,llm,1,1) - call resume_timer(timer_dissip) -c$OMP END MASTER -c$OMP BARRIER - CALL grad_p (klevel,divgra, ghx , ghy ) - -c$OMP BARRIER -c$OMP MASTER - call suspend_timer(timer_dissip) - call exchange_Hallo(ghy,ip1jm,llm,1,0) - call resume_timer(timer_dissip) -c$OMP END MASTER -c$OMP BARRIER - - CALL diverg_p (klevel, ghx , ghy , divgra ) - - jjb=jj_begin - jje=jj_end - CALL filtreg_p( divgra,jjb,jje,jjp1,klevel,2,1,.true.,1) - -c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) - DO 5 l = 1,klevel - DO 4 ij = ijb, ije - divgra_out( ij,l ) = - cdivh * divgra( ij,l ) - 4 CONTINUE - 5 CONTINUE -c$OMP END DO NOWAIT -c - 10 CONTINUE - RETURN - END Index: LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc.F (revision 1748) +++ LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc.F (revision 1749) @@ -27,4 +27,5 @@ c USE parallel + USE mod_filtreg_p IMPLICIT NONE c @@ -120,5 +121,6 @@ jjb=jj_begin jje=jj_end - CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 ) + CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, + & 2, 1, .TRUE., 1 ) ! our work is done, exit routine @@ -206,5 +208,6 @@ jjb=jj_begin jje=jj_end - CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 ) + CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, + & 2, 1, .TRUE., 1 ) c EST-CE UTILE ?? : calcul de beta Index: LMDZ5/trunk/libf/dyn3dmem/gcm.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/gcm.F (revision 1748) +++ LMDZ5/trunk/libf/dyn3dmem/gcm.F (revision 1749) @@ -270,5 +270,5 @@ ! constants & fields, if we run the 'newtonian' or 'SW' cases: if (iflag_phys.ne.1) then - CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) + CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0) endif @@ -291,5 +291,5 @@ . 'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT' if (.not.read_start) then - CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) + CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0) endif @@ -398,5 +398,7 @@ #endif - + if (iflag_phys.eq.1) then + ! these initialisations have already been done (via iniacademic) + ! if running in SW or Newtonian mode c----------------------------------------------------------------------- c Initialisation des constantes dynamiques : @@ -414,4 +416,5 @@ c -------------------------- CALL inifilr + endif ! of if (iflag_phys.eq.1) c c----------------------------------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dmem/gradiv_p.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/gradiv_p.F (revision 1748) +++ (revision ) @@ -1,109 +1,0 @@ - SUBROUTINE gradiv_p(klevel, xcov, ycov, ld, gdx_out, gdy_out ) -c -c Auteur : P. Le Van -c -c *************************************************************** -c -c ld -c calcul de (grad (div) ) du vect. v .... -c -c xcov et ycov etant les composant.covariantes de v -c **************************************************************** -c xcov , ycov et ld sont des arguments d'entree pour le s-prog -c gdx et gdy sont des arguments de sortie pour le s-prog -c -c - USE parallel - USE times - IMPLICIT NONE -c -#include "dimensions.h" -#include "paramet.h" -#include "comdissipn.h" -#include "logic.h" - - INTEGER klevel -c - REAL xcov( ip1jmp1,klevel ), ycov( ip1jm,klevel ) - REAL,SAVE :: gdx( ip1jmp1,llm ), gdy( ip1jm,llm ) - - REAL gdx_out( ip1jmp1,klevel ), gdy_out( ip1jm,klevel ) - - REAL,SAVE :: div(ip1jmp1,llm) - - INTEGER l,ij,iter,ld -c - INTEGER ijb,ije,jjb,jje -c -c -c CALL SCOPY( ip1jmp1*klevel,xcov,1,gdx,1 ) -c CALL SCOPY( ip1jm*klevel, ycov,1,gdy,1 ) - - ijb=ij_begin - ije=ij_end - -c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) - DO l = 1,klevel - gdx(ijb:ije,l)=xcov(ijb:ije,l) - ENDDO -c$OMP END DO NOWAIT - - ijb=ij_begin - ije=ij_end - if(pole_sud) ije=ij_end-iip1 - -c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) - DO l = 1,klevel - gdy(ijb:ije,l)=ycov(ijb:ije,l) - ENDDO -c$OMP END DO NOWAIT - -c - DO 10 iter = 1,ld - -c$OMP BARRIER -c$OMP MASTER - call suspend_timer(timer_dissip) - call exchange_Hallo(gdy,ip1jm,llm,1,0) - call resume_timer(timer_dissip) -c$OMP END MASTER -c$OMP BARRIER - - CALL diverg_p( klevel, gdx , gdy, div ) - - jjb=jj_begin - jje=jj_end - CALL filtreg_p( div,jjb,jje, jjp1, klevel, 2,1, .true.,2 ) - -c call exchange_Hallo(div,ip1jmp1,llm,0,1) - -c$OMP BARRIER -c$OMP MASTER - call suspend_timer(timer_dissip) - call exchange_Hallo(div,ip1jmp1,llm,1,1) - call resume_timer(timer_dissip) -c$OMP END MASTER -c$OMP BARRIER - - CALL grad_p( klevel, div, gdx, gdy ) -c - -c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) - DO 5 l = 1, klevel - - if(pole_sud) ije=ij_end - DO 3 ij = ijb, ije - gdx_out( ij,l ) = - gdx( ij,l ) * cdivu - 3 CONTINUE - - if(pole_sud) ije=ij_end-iip1 - DO 4 ij = ijb, ije - gdy_out( ij,l ) = - gdy( ij,l ) * cdivu - 4 CONTINUE - - 5 CONTINUE -c$OMP END DO NOWAIT -c - 10 CONTINUE - RETURN - END Index: LMDZ5/trunk/libf/dyn3dmem/iniacademic.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/iniacademic.F90 (revision 1748) +++ (revision ) @@ -1,277 +1,0 @@ -! -! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $ -! -SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0) - - USE filtreg_mod - USE infotrac, ONLY : nqtot - USE control_mod, ONLY: day_step,planet_type -#ifdef CPP_IOIPSL - USE IOIPSL -#else - ! if not using IOIPSL, we still need to use (a local version of) getin - USE ioipsl_getincom -#endif - USE Write_Field - - ! Author: Frederic Hourdin original: 15/01/93 - ! The forcing defined here is from Held and Suarez, 1994, Bulletin - ! of the American Meteorological Society, 75, 1825. - - IMPLICIT NONE - - ! Declararations: - ! --------------- - - include "dimensions.h" - include "paramet.h" - include "comvert.h" - include "comconst.h" - include "comgeom.h" - include "academic.h" - include "ener.h" - include "temps.h" - include "iniprint.h" - include "logic.h" - - ! Arguments: - ! ---------- - - real time_0 - - ! variables dynamiques - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants - REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL q(ip1jmp1,llm,nqtot) ! champs advectes - REAL ps(ip1jmp1) ! pression au sol - REAL masse(ip1jmp1,llm) ! masse d'air - REAL phis(ip1jmp1) ! geopotentiel au sol - - ! Local: - ! ------ - - REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches - REAL pks(ip1jmp1) ! exner au sol - REAL pk(ip1jmp1,llm) ! exner au milieu des couches - REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches - REAL phi(ip1jmp1,llm) ! geopotentiel - REAL ddsin,zsig,tetapv,w_pv ! variables auxiliaires - real tetastrat ! potential temperature in the stratosphere, in K - real tetajl(jjp1,llm) - INTEGER i,j,l,lsup,ij - - REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T - REAL k_f,k_c_a,k_c_s ! Constantes de rappel - LOGICAL ok_geost ! Initialisation vent geost. ou nul - LOGICAL ok_pv ! Polar Vortex - REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex - - real zz,ran1 - integer idum - - REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr - - character(len=*),parameter :: modname="iniacademic" - character(len=80) :: abort_message - - !----------------------------------------------------------------------- - ! 1. Initializations for Earth-like case - ! -------------------------------------- - ! - ! initialize planet radius, rotation rate,... - call conf_planete - - time_0=0. - day_ref=1 - annee_ref=0 - - im = iim - jm = jjm - day_ini = 1 - dtvr = daysec/REAL(day_step) - zdtvr=dtvr - etot0 = 0. - ptot0 = 0. - ztot0 = 0. - stot0 = 0. - ang0 = 0. - - if (llm == 1) then - ! specific initializations for the shallow water case - kappa=1 - endif - - CALL iniconst - CALL inigeom - CALL inifilr - - if (llm == 1) then - ! initialize fields for the shallow water case, if required - if (.not.read_start) then - phis(:)=0. - q(:,:,:)=0 - CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps) - endif - endif - - academic_case: if (iflag_phys >= 2) then - ! initializations - - ! 1. local parameters - ! by convention, winter is in the southern hemisphere - ! Geostrophic wind or no wind? - ok_geost=.TRUE. - CALL getin('ok_geost',ok_geost) - ! Constants for Newtonian relaxation and friction - k_f=1. !friction - CALL getin('k_j',k_f) - k_f=1./(daysec*k_f) - k_c_s=4. !cooling surface - CALL getin('k_c_s',k_c_s) - k_c_s=1./(daysec*k_c_s) - k_c_a=40. !cooling free atm - CALL getin('k_c_a',k_c_a) - k_c_a=1./(daysec*k_c_a) - ! Constants for Teta equilibrium profile - teta0=315. ! mean Teta (S.H. 315K) - CALL getin('teta0',teta0) - ttp=200. ! Tropopause temperature (S.H. 200K) - CALL getin('ttp',ttp) - eps=0. ! Deviation to N-S symmetry(~0-20K) - CALL getin('eps',eps) - delt_y=60. ! Merid Temp. Gradient (S.H. 60K) - CALL getin('delt_y',delt_y) - delt_z=10. ! Vertical Gradient (S.H. 10K) - CALL getin('delt_z',delt_z) - ! Polar vortex - ok_pv=.false. - CALL getin('ok_pv',ok_pv) - phi_pv=-50. ! Latitude of edge of vortex - CALL getin('phi_pv',phi_pv) - phi_pv=phi_pv*pi/180. - dphi_pv=5. ! Width of the edge - CALL getin('dphi_pv',dphi_pv) - dphi_pv=dphi_pv*pi/180. - gam_pv=4. ! -dT/dz vortex (in K/km) - CALL getin('gam_pv',gam_pv) - - ! 2. Initialize fields towards which to relax - ! Friction - knewt_g=k_c_a - DO l=1,llm - zsig=presnivs(l)/preff - knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) - kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) - ENDDO - DO j=1,jjp1 - clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 - ENDDO - - ! Potential temperature - DO l=1,llm - zsig=presnivs(l)/preff - tetastrat=ttp*zsig**(-kappa) - tetapv=tetastrat - IF ((ok_pv).AND.(zsig.LT.0.1)) THEN - tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) - ENDIF - DO j=1,jjp1 - ! Troposphere - ddsin=sin(rlatu(j)) - tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & - -delt_z*(1.-ddsin*ddsin)*log(zsig) - if (planet_type=="giant") then - tetajl(j,l)=teta0+(delt_y* & - ((sin(rlatu(j)*3.14159*eps+0.0001))**2) & - / ((rlatu(j)*3.14159*eps+0.0001)**2)) & - -delt_z*log(zsig) - endif - ! Profil stratospherique isotherme (+vortex) - w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. - tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv - tetajl(j,l)=MAX(tetajl(j,l),tetastrat) - ENDDO - ENDDO - - ! CALL writefield('theta_eq',tetajl) - - do l=1,llm - do j=1,jjp1 - do i=1,iip1 - ij=(j-1)*iip1+i - tetarappel(ij,l)=tetajl(j,l) - enddo - enddo - enddo - - ! 3. Initialize fields (if necessary) - IF (.NOT. read_start) THEN - ! surface pressure - if (iflag_phys>2) then - ! specific value for CMIP5 aqua/terra planets - ! "Specify the initial dry mass to be equivalent to - ! a global mean surface pressure (101325 minus 245) Pa." - ps(:)=101080. - else - ! use reference surface pressure - ps(:)=preff - endif - - ! ground geopotential - phis(:)=0. - - CALL pression ( ip1jmp1, ap, bp, ps, p ) - if (pressure_exner) then - CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf ) - else - call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf) - endif - CALL massdair(p,masse) - - ! bulk initialization of temperature - teta(:,:)=tetarappel(:,:) - - ! geopotential - CALL geopot(ip1jmp1,teta,pk,pks,phis,phi) - - ! winds - if (ok_geost) then - call ugeostr(phi,ucov) - else - ucov(:,:)=0. - endif - vcov(:,:)=0. - - ! bulk initialization of tracers - if (planet_type=="earth") then - ! Earth: first two tracers will be water - do i=1,nqtot - if (i == 1) q(:,:,i)=1.e-10 - if (i == 2) q(:,:,i)=1.e-15 - if (i.gt.2) q(:,:,i)=0. - enddo - else - q(:,:,:)=0 - endif ! of if (planet_type=="earth") - - ! add random perturbation to temperature - idum = -1 - zz = ran1(idum) - idum = 0 - do l=1,llm - do ij=iip2,ip1jm - teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum)) - enddo - enddo - - ! maintain periodicity in longitude - do l=1,llm - do ij=1,ip1jmp1,iip1 - teta(ij+iim,l)=teta(ij,l) - enddo - enddo - - ENDIF ! of IF (.NOT. read_start) - endif academic_case - -END SUBROUTINE iniacademic Index: LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90 =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90 (revision 1749) +++ LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90 (revision 1749) @@ -0,0 +1,308 @@ +! +! $Id: iniacademic.F90 1625 2012-05-09 13:14:48Z lguez $ +! +SUBROUTINE iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0) + + USE filtreg_mod + USE infotrac, ONLY : nqtot + USE control_mod, ONLY: day_step,planet_type + USE parallel +#ifdef CPP_IOIPSL + USE IOIPSL +#else + ! if not using IOIPSL, we still need to use (a local version of) getin + USE ioipsl_getincom +#endif + USE Write_Field + + ! Author: Frederic Hourdin original: 15/01/93 + ! The forcing defined here is from Held and Suarez, 1994, Bulletin + ! of the American Meteorological Society, 75, 1825. + + IMPLICIT NONE + + ! Declararations: + ! --------------- + + include "dimensions.h" + include "paramet.h" + include "comvert.h" + include "comconst.h" + include "comgeom.h" + include "academic.h" + include "ener.h" + include "temps.h" + include "iniprint.h" + include "logic.h" + + ! Arguments: + ! ---------- + + real time_0 + + ! variables dynamiques + REAL vcov(ijb_v:ije_v,llm),ucov(ijb_u:ije_u,llm) ! vents covariants + REAL teta(ijb_u:ije_u,llm) ! temperature potentielle + REAL q(ijb_u:ije_u,llm,nqtot) ! champs advectes + REAL ps(ijb_u:ije_u) ! pression au sol + REAL masse(ijb_u:ije_u,llm) ! masse d'air + REAL phis(ijb_u:ije_u) ! geopotentiel au sol + + ! Local: + ! ------ + + REAL,ALLOCATABLE :: vcov_glo(:,:),ucov_glo(:,:),teta_glo(:,:) + REAL,ALLOCATABLE :: q_glo(:,:),masse_glo(:,:),ps_glo(:) + REAL,ALLOCATABLE :: phis_glo(:) + REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches + REAL pks(ip1jmp1) ! exner au sol + REAL pk(ip1jmp1,llm) ! exner au milieu des couches + REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches + REAL phi(ip1jmp1,llm) ! geopotentiel + REAL ddsin,zsig,tetapv,w_pv ! variables auxiliaires + real tetastrat ! potential temperature in the stratosphere, in K + real tetajl(jjp1,llm) + INTEGER i,j,l,lsup,ij + + REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T + REAL k_f,k_c_a,k_c_s ! Constantes de rappel + LOGICAL ok_geost ! Initialisation vent geost. ou nul + LOGICAL ok_pv ! Polar Vortex + REAL phi_pv,dphi_pv,gam_pv ! Constantes pour polar vortex + + real zz,ran1 + integer idum + + REAL zdtvr + real,allocatable :: alpha(:,:),beta(:,:) + + character(len=*),parameter :: modname="iniacademic" + character(len=80) :: abort_message + + !----------------------------------------------------------------------- + ! 1. Initializations for Earth-like case + ! -------------------------------------- + ! + ! initialize planet radius, rotation rate,... + call conf_planete + + time_0=0. + day_ref=1 + annee_ref=0 + + im = iim + jm = jjm + day_ini = 1 + dtvr = daysec/REAL(day_step) + zdtvr=dtvr + etot0 = 0. + ptot0 = 0. + ztot0 = 0. + stot0 = 0. + ang0 = 0. + + if (llm == 1) then + ! specific initializations for the shallow water case + kappa=1 + endif + + CALL iniconst + CALL inigeom + CALL inifilr + + if (llm == 1) then + ! initialize fields for the shallow water case, if required + if (.not.read_start) then + phis(ijb_u:ije_u)=0. + q(ijb_u:ije_u,1:llm,1:nqtot)=0 + CALL sw_case_williamson91_6_loc(vcov,ucov,teta,masse,ps) + endif + endif + + academic_case: if (iflag_phys >= 2) then + ! initializations + + ! 1. local parameters + ! by convention, winter is in the southern hemisphere + ! Geostrophic wind or no wind? + ok_geost=.TRUE. + CALL getin('ok_geost',ok_geost) + ! Constants for Newtonian relaxation and friction + k_f=1. !friction + CALL getin('k_j',k_f) + k_f=1./(daysec*k_f) + k_c_s=4. !cooling surface + CALL getin('k_c_s',k_c_s) + k_c_s=1./(daysec*k_c_s) + k_c_a=40. !cooling free atm + CALL getin('k_c_a',k_c_a) + k_c_a=1./(daysec*k_c_a) + ! Constants for Teta equilibrium profile + teta0=315. ! mean Teta (S.H. 315K) + CALL getin('teta0',teta0) + ttp=200. ! Tropopause temperature (S.H. 200K) + CALL getin('ttp',ttp) + eps=0. ! Deviation to N-S symmetry(~0-20K) + CALL getin('eps',eps) + delt_y=60. ! Merid Temp. Gradient (S.H. 60K) + CALL getin('delt_y',delt_y) + delt_z=10. ! Vertical Gradient (S.H. 10K) + CALL getin('delt_z',delt_z) + ! Polar vortex + ok_pv=.false. + CALL getin('ok_pv',ok_pv) + phi_pv=-50. ! Latitude of edge of vortex + CALL getin('phi_pv',phi_pv) + phi_pv=phi_pv*pi/180. + dphi_pv=5. ! Width of the edge + CALL getin('dphi_pv',dphi_pv) + dphi_pv=dphi_pv*pi/180. + gam_pv=4. ! -dT/dz vortex (in K/km) + CALL getin('gam_pv',gam_pv) + + ! 2. Initialize fields towards which to relax + ! Friction + knewt_g=k_c_a + DO l=1,llm + zsig=presnivs(l)/preff + knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3) + kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3) + ENDDO + DO j=1,jjp1 + clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4 + ENDDO + + ! Potential temperature + DO l=1,llm + zsig=presnivs(l)/preff + tetastrat=ttp*zsig**(-kappa) + tetapv=tetastrat + IF ((ok_pv).AND.(zsig.LT.0.1)) THEN + tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g) + ENDIF + DO j=1,jjp1 + ! Troposphere + ddsin=sin(rlatu(j)) + tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin & + -delt_z*(1.-ddsin*ddsin)*log(zsig) + if (planet_type=="giant") then + tetajl(j,l)=teta0+(delt_y* & + ((sin(rlatu(j)*3.14159*eps+0.0001))**2) & + / ((rlatu(j)*3.14159*eps+0.0001)**2)) & + -delt_z*log(zsig) + endif + ! Profil stratospherique isotherme (+vortex) + w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2. + tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv + tetajl(j,l)=MAX(tetajl(j,l),tetastrat) + ENDDO + ENDDO + + ! CALL writefield('theta_eq',tetajl) + + do l=1,llm + do j=1,jjp1 + do i=1,iip1 + ij=(j-1)*iip1+i + tetarappel(ij,l)=tetajl(j,l) + enddo + enddo + enddo + + ! 3. Initialize fields (if necessary) + IF (.NOT. read_start) THEN + ! allocate global fields: +! allocate(vcov_glo(ip1jm,llm)) + allocate(ucov_glo(ip1jmp1,llm)) + allocate(teta_glo(ip1jmp1,llm)) + allocate(ps_glo(ip1jmp1)) + allocate(masse_glo(ip1jmp1,llm)) + allocate(phis_glo(ip1jmp1)) + allocate(alpha(ip1jmp1,llm)) + allocate(beta(ip1jmp1,llm)) + + ! surface pressure + if (iflag_phys>2) then + ! specific value for CMIP5 aqua/terra planets + ! "Specify the initial dry mass to be equivalent to + ! a global mean surface pressure (101325 minus 245) Pa." + ps_glo(:)=101080. + else + ! use reference surface pressure + ps_glo(:)=preff + endif + + ! ground geopotential + phis_glo(:)=0. + + CALL pression ( ip1jmp1, ap, bp, ps_glo, p ) + if (pressure_exner) then + CALL exner_hyb( ip1jmp1, ps_glo, p,alpha,beta, pks, pk, pkf ) + else + call exner_milieu(ip1jmp1,ps_glo,p,beta,pks,pk,pkf) + endif + CALL massdair(p,masse_glo) + + ! bulk initialization of temperature + teta_glo(:,:)=tetarappel(:,:) + + ! geopotential + CALL geopot(ip1jmp1,teta_glo,pk,pks,phis_glo,phi) + + ! winds + if (ok_geost) then + call ugeostr(phi,ucov_glo) + else + ucov_glo(:,:)=0. + endif + vcov(ijb_v:ije_v,1:llm)=0. + + ! bulk initialization of tracers + if (planet_type=="earth") then + ! Earth: first two tracers will be water + do i=1,nqtot + if (i == 1) q(ijb_u:ije_u,:,i)=1.e-10 + if (i == 2) q(ijb_u:ije_u,:,i)=1.e-15 + if (i.gt.2) q(ijb_u:ije_u,:,i)=0. + enddo + else + q(ijb_u:ije_u,:,:)=0 + endif ! of if (planet_type=="earth") + + ! add random perturbation to temperature + idum = -1 + zz = ran1(idum) + idum = 0 + do l=1,llm + do ij=iip2,ip1jm + teta_glo(ij,l)=teta_glo(ij,l)*(1.+0.005*ran1(idum)) + enddo + enddo + + ! maintain periodicity in longitude + do l=1,llm + do ij=1,ip1jmp1,iip1 + teta_glo(ij+iim,l)=teta_glo(ij,l) + enddo + enddo + + ! copy data from global array to local array: + teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:) + ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:) +! vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:) + masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:) + ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u) + phis(ijb_u:ije_u)=phis_glo(ijb_u:ije_u) + + deallocate(teta_glo) + deallocate(ucov_glo) +! deallocate(vcov_glo) + deallocate(masse_glo) + deallocate(ps_glo) + deallocate(phis_glo) + deallocate(alpha) + deallocate(beta) + ENDIF ! of IF (.NOT. read_start) + endif academic_case + +END SUBROUTINE iniacademic_loc Index: LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6.F (revision 1748) +++ (revision ) @@ -1,140 +1,0 @@ -! -! $Id $ -! - SUBROUTINE sw_case_williamson91_6(vcov,ucov,teta,masse,ps) - -c======================================================================= -c -c Author: Thomas Dubos original: 26/01/2010 -c ------- -c -c Subject: -c ------ -c Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz -c -c Method: -c -------- -c -c Interface: -c ---------- -c -c Input: -c ------ -c -c Output: -c ------- -c -c======================================================================= - IMPLICIT NONE -c----------------------------------------------------------------------- -c Declararations: -c --------------- - -#include "dimensions.h" -#include "paramet.h" -#include "comvert.h" -#include "comconst.h" -#include "comgeom.h" -#include "iniprint.h" - -c Arguments: -c ---------- - -c variables dynamiques - REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants - REAL teta(ip1jmp1,llm) ! temperature potentielle - REAL ps(ip1jmp1) ! pression au sol - REAL masse(ip1jmp1,llm) ! masse d'air - REAL phis(ip1jmp1) ! geopotentiel au sol - -c Local: -c ------ - - REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches - REAL pks(ip1jmp1) ! exner au sol - REAL pk(ip1jmp1,llm) ! exner au milieu des couches - REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches - REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) - - REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps - INTEGER i,j,ij - - REAL, PARAMETER :: rho=1 ! masse volumique de l'air (arbitraire) - REAL, PARAMETER :: K = 7.848e-6 ! K = \omega - REAL, PARAMETER :: gh0 = 9.80616 * 8e3 - INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2 ! mode 4 -c NB : rad = 6371220 dans W91 (6371229 dans LMDZ) -c omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ) - - IF(0==0) THEN -c Williamson et al. (1991) : onde de Rossby-Haurwitz - teta = preff/rho/cpp -c geopotentiel (pression de surface) - do j=1,jjp1 - costh2 = cos(rlatu(j))**2 - Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0 - Ath = .25*(K**2)*(costh2**(R0-1))*Ath - Ath = .5*K*(2*omeg+K)*costh2 + Ath - Bth = (R1*R1+1)-R1*R1*costh2 - Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth - Cth = R1*costh2 - R2 - Cth = .25*K*K*(costh2**R0)*Cth - do i=1,iip1 - ij=(j-1)*iip1+i - lon = rlonv(i) - dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon) - ps(ij) = rho*(gh0 + (rad**2)*dps) - enddo - enddo - write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps) -c vitesse zonale ucov - do j=1,jjp1 - costh = cos(rlatu(j)) - costh2 = costh**2 - Ath = rad*K*costh - Bth = R0*(1-costh2)-costh2 - Bth = rad*K*Bth*(costh**(R0-1)) - do i=1,iip1 - ij=(j-1)*iip1+i - lon = rlonu(i) - ucov(ij,1) = (Ath + Bth*cos(R0*lon)) - enddo - enddo - write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1)) - ucov(:,1)=ucov(:,1)*cu -c vitesse meridienne vcov - do j=1,jjm - sinth = sin(rlatv(j)) - costh = cos(rlatv(j)) - Ath = -rad*K*R0*sinth*(costh**(R0-1)) - do i=1,iip1 - ij=(j-1)*iip1+i - lon = rlonv(i) - vcov(ij,1) = Ath*sin(R0*lon) - enddo - enddo - write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1)) - vcov(:,1)=vcov(:,1)*cv - -c ucov=0 -c vcov=0 - ELSE -c test non-tournant, onde se propageant en latitude - do j=1,jjp1 - do i=1,iip1 - ij=(j-1)*iip1+i - ps(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2) ) - enddo - enddo - -c rho = preff/(cpp*teta) - teta = .01*preff/cpp ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j - ucov=0. - vcov=0. - END IF - - CALL pression ( ip1jmp1, ap, bp, ps, p ) - CALL massdair(p,masse) - - END -c----------------------------------------------------------------------- Index: LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6_loc.F =================================================================== --- LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6_loc.F (revision 1749) +++ LMDZ5/trunk/libf/dyn3dmem/sw_case_williamson91_6_loc.F (revision 1749) @@ -0,0 +1,189 @@ +! +! $Id $ +! + SUBROUTINE sw_case_williamson91_6_loc(vcov,ucov,teta,masse,ps) + +c======================================================================= +c +c Author: Thomas Dubos original: 26/01/2010 +c ------- +c +c Subject: +c ------ +c Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz +c +c Method: +c -------- +c +c Interface: +c ---------- +c +c Input: +c ------ +c +c Output: +c ------- +c +c======================================================================= + USE parallel + + IMPLICIT NONE +c----------------------------------------------------------------------- +c Declararations: +c --------------- + +#include "dimensions.h" +#include "paramet.h" +#include "comvert.h" +#include "comconst.h" +#include "comgeom.h" +#include "iniprint.h" + +c Arguments: +c ---------- + +c variables dynamiques + REAL vcov(ijb_v:ije_v,llm),ucov(ijb_u:ije_u,llm) ! vents covariants + REAL teta(ijb_u:ije_u,llm) ! temperature potentielle + REAL ps(ijb_u:ije_u) ! pression au sol + REAL masse(ijb_u:ije_u,llm) ! masse d'air + REAL phis(ijb_u:ije_u) ! geopotentiel au sol + +c Local: +c ------ + + real,allocatable :: ucov_glo(:,:) + real,allocatable :: vcov_glo(:,:) + real,allocatable :: teta_glo(:,:) + real,allocatable :: masse_glo(:,:) + real,allocatable :: ps_glo(:) + +! REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches +! REAL pks(ip1jmp1) ! exner au sol +! REAL pk(ip1jmp1,llm) ! exner au milieu des couches +! REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches +! REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm) + + real,allocatable :: p(:,:) + real,allocatable :: pks(:) + real,allocatable :: pk(:,:) + real,allocatable :: pkf(:,:) + real,allocatable :: alpha(:,:),beta(:,:) + + REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps + INTEGER i,j,ij + + REAL, PARAMETER :: rho=1 ! masse volumique de l'air (arbitraire) + REAL, PARAMETER :: K = 7.848e-6 ! K = \omega + REAL, PARAMETER :: gh0 = 9.80616 * 8e3 + INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2 ! mode 4 +c NB : rad = 6371220 dans W91 (6371229 dans LMDZ) +c omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ) + + + ! allocate (global) arrays + allocate(vcov_glo(ip1jm,llm)) + allocate(ucov_glo(ip1jmp1,llm)) + allocate(teta_glo(ip1jmp1,llm)) + allocate(ps_glo(ip1jmp1)) + allocate(masse_glo(ip1jmp1,llm)) + + allocate(p(ip1jmp1,llmp1)) + allocate(pks(ip1jmp1)) + allocate(pk(ip1jmp1,llm)) + allocate(pkf(ip1jmp1,llm)) + allocate(alpha(ip1jmp1,llm)) + allocate(beta(ip1jmp1,llm)) + + IF(0==0) THEN +!c Williamson et al. (1991) : onde de Rossby-Haurwitz + teta_glo(:,:) = preff/rho/cpp +!c geopotentiel (pression de surface) + do j=1,jjp1 + costh2 = cos(rlatu(j))**2 + Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0 + Ath = .25*(K**2)*(costh2**(R0-1))*Ath + Ath = .5*K*(2*omeg+K)*costh2 + Ath + Bth = (R1*R1+1)-R1*R1*costh2 + Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth + Cth = R1*costh2 - R2 + Cth = .25*K*K*(costh2**R0)*Cth + do i=1,iip1 + ij=(j-1)*iip1+i + lon = rlonv(i) + dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon) + ps_glo(ij) = rho*(gh0 + (rad**2)*dps) + enddo + enddo +! write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps) +c vitesse zonale ucov + do j=1,jjp1 + costh = cos(rlatu(j)) + costh2 = costh**2 + Ath = rad*K*costh + Bth = R0*(1-costh2)-costh2 + Bth = rad*K*Bth*(costh**(R0-1)) + do i=1,iip1 + ij=(j-1)*iip1+i + lon = rlonu(i) + ucov_glo(ij,1) = (Ath + Bth*cos(R0*lon)) + enddo + enddo +! write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1)) + ucov_glo(:,1)=ucov_glo(:,1)*cu +c vitesse meridienne vcov + do j=1,jjm + sinth = sin(rlatv(j)) + costh = cos(rlatv(j)) + Ath = -rad*K*R0*sinth*(costh**(R0-1)) + do i=1,iip1 + ij=(j-1)*iip1+i + lon = rlonv(i) + vcov_glo(ij,1) = Ath*sin(R0*lon) + enddo + enddo + write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1)) + vcov_glo(:,1)=vcov_glo(:,1)*cv + +c ucov_glo=0 +c vcov_glo=0 + ELSE +c test non-tournant, onde se propageant en latitude + do j=1,jjp1 + do i=1,iip1 + ij=(j-1)*iip1+i + ps_glo(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2)) + enddo + enddo + +c rho = preff/(cpp*teta) + teta_glo(:,:) = .01*preff/cpp ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j + ucov_glo(:,:)=0. + vcov_glo(:,:)=0. + END IF + + CALL pression ( ip1jmp1, ap, bp, ps_glo, p ) + CALL massdair(p,masse_glo) + + ! copy data from global array to local array: + teta(ijb_u:ije_u,:)=teta_glo(ijb_u:ije_u,:) + ucov(ijb_u:ije_u,:)=ucov_glo(ijb_u:ije_u,:) + vcov(ijb_v:ije_v,:)=vcov_glo(ijb_v:ije_v,:) + masse(ijb_u:ije_u,:)=masse_glo(ijb_u:ije_u,:) + ps(ijb_u:ije_u)=ps_glo(ijb_u:ije_u) + + ! cleanup + deallocate(teta_glo) + deallocate(ucov_glo) + deallocate(vcov_glo) + deallocate(masse_glo) + deallocate(ps_glo) + deallocate(p) + deallocate(pks) + deallocate(pk) + deallocate(pkf) + deallocate(alpha) + deallocate(beta) + + END +c-----------------------------------------------------------------------