MODULE etat0dyn ! !******************************************************************************* ! Purpose: Create dynamical initial state using atmospheric fields from a ! database of atmospheric to initialize the model. !------------------------------------------------------------------------------- ! Comments: ! ! * This module is designed to work for Earth (and with ioipsl) ! ! * etat0dyn_netcdf routine can access to NetCDF data through the following ! routine (to be called after restget): ! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,& ! champ, val_exp, lon_in2, lat_in2, ibar) ! ! * Variables should have the following names in the NetCDF files: ! 'U' : East ward wind (in "ECDYN.nc") ! 'V' : Northward wind (in "ECDYN.nc") ! 'TEMP' : Temperature (in "ECDYN.nc") ! 'R' : Relative humidity (in "ECDYN.nc") ! 'RELIEF' : High resolution orography (in "Relief.nc") ! ! * The land mask and corresponding weights can be: ! 1) already known (in particular if etat0dyn has been called before) ; ! in this case, ANY(masque(:,:)/=-99999.) = .TRUE. ! 2) computed using the ocean mask from the ocean model (to ensure ocean ! fractions are the same for atmosphere and ocean) for coupled runs. ! File name: "o2a.nc" ; variable name: "OceMask" ! 3) computed from topography file "Relief.nc" for forced runs. ! ! * There is a big mess with the longitude size. Should it be iml or iml+1 ? ! I have chosen to use the iml+1 as an argument to this routine and we declare ! internaly smaller fields when needed. This needs to be cleared once and for ! all in LMDZ. A convention is required. !------------------------------------------------------------------------------- USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo USE assert_eq_m, ONLY: assert_eq #ifdef CPP_PHYS USE indice_sol_mod, ONLY: epsfra #endif IMPLICIT NONE PRIVATE PUBLIC :: etat0dyn_netcdf include "iniprint.h" include "dimensions.h" include "paramet.h" include "comgeom2.h" include "comvert.h" include "comconst.h" include "temps.h" include "comdissnew.h" include "serre.h" REAL, SAVE :: deg2rad #ifndef CPP_PHYS REAL, SAVE :: epsfra= 1.E-5 #endif INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:) CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc' CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc' CONTAINS !------------------------------------------------------------------------------- ! SUBROUTINE etat0dyn_netcdf(ib, masque, phis) ! !------------------------------------------------------------------------------- ! Purpose: Create dynamical initial states. !------------------------------------------------------------------------------- ! Notes: 1) This routine is designed to work for Earth ! 2) If masque(:,:)/=-99999., masque and phis are already known. ! Otherwise: read it from ocean model file (coupled run) or compute it. !------------------------------------------------------------------------------- USE control_mod #ifdef CPP_PHYS #ifdef CPP_EARTH USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz USE regr_pr_o3_m, ONLY: regr_pr_o3 USE press_coefoz_m, ONLY: press_coefoz #endif #endif USE exner_hyb_m, ONLY: exner_hyb USE exner_milieu_m, ONLY: exner_milieu USE infotrac USE filtreg_mod !#endif IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: LOGICAL, INTENT(IN) :: ib !--- Barycentric interpolation REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential !#ifndef CPP_EARTH ! WRITE(lunout,*)'etat0dyn_netcdf: Earth-specific routine, needs Earth physics' !#else !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: modname, fmt INTEGER :: i, j, l, ji, itau, iday REAL :: xpn, xps, time, phystep REAL, DIMENSION(iip1,jjp1) :: psol, masque_tmp REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd REAL, DIMENSION(iip1,jjp1,llm) :: pk, pls, y, masse REAL, DIMENSION(iip1,jjm ,llm) :: vvent REAL, DIMENSION(ip1jm ,llm) :: pbarv REAL, DIMENSION(ip1jmp1 ,llm) :: pbaru, phi, w REAL, DIMENSION(ip1jmp1) :: pks REAL, DIMENSION(iim) :: xppn, xpps REAL, ALLOCATABLE :: q3d(:,:,:,:) !------------------------------------------------------------------------------- modname='etat0dyn_netcdf' deg2rad = pi/180.0 ! Initializations for tracers and filter !******************************************************************************* CALL infotrac_init ALLOCATE(q3d(iip1,jjp1,llm,nqtot)) CALL inifilr() ! Compute ground geopotential and possibly the mask. !******************************************************************************* masque_tmp(:,:)=masque(:,:) CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp) WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt) IF(ALL(masque==-99999.)) THEN !--- KEEP NEW MASK masque=masque_tmp IF(prt_level>=1) THEN WRITE(lunout,*)'BUILT MASK :' WRITE(lunout,fmt) NINT(masque) END IF WHERE( masque(:,:)mass fraction END IF #endif ! Intermediate computation !******************************************************************************* CALL massdair(p3d,masse) WRITE(lunout,*)' ALPHAX ',alphax DO l=1,llm xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l) xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l) xpn=SUM(xppn)/apoln xps=SUM(xpps)/apols masse(:,1 ,l)=xpn masse(:,jjp1,l)=xps END DO q3d(iip1,:,:,:)=q3d(1,:,:,:) ! Writing !******************************************************************************* CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & tetatemp, vert_prof_dissip) WRITE(lunout,*)'sortie inidissip' itau=0 itau_dyn=0 itau_phy=0 iday=dayref+itau/day_step time=FLOAT(itau-(iday-dayref)*day_step)/day_step IF(time>1.) THEN time=time-1 iday=iday+1 END IF day_ref=dayref annee_ref=anneeref CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) WRITE(lunout,*)'sortie geopot' CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & phi, w, pbaru, pbarv, time+iday-dayref) WRITE(lunout,*)'sortie caldyn0' CALL dynredem0( "start.nc", dayref, phis) WRITE(lunout,*)'sortie dynredem0' CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) WRITE(lunout,*)'sortie dynredem1' CALL histclo() !#endif !#endif of #ifdef CPP_EARTH END SUBROUTINE etat0dyn_netcdf ! !------------------------------------------------------------------------------- !#ifdef CPP_EARTH !------------------------------------------------------------------------------- ! SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,& champ, val_exp, lon_in2, lat_in2, ibar) !------------------------------------------------------------------------------- IMPLICIT NONE !=============================================================================== ! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R ! (3D fields) of file dynfname. !------------------------------------------------------------------------------- ! Note: An input auxilliary field "workvar" has to be specified in two cases: ! * for "q": the saturated humidity. ! * for "topot": the Exner function. !=============================================================================== ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: var REAL, INTENT(IN) :: lon_in(:) ! dim (iml) REAL, INTENT(IN) :: lat_in(:) ! dim (jml) REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml) REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml) REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml) REAL, INTENT(IN) :: val_exp REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) LOGICAL, INTENT(IN) :: ibar !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=10) :: vname CHARACTER(LEN=256) :: msg, modname="startget_dyn3d" INTEGER :: iml, jml, jml2, lml, il REAL :: xppn, xpps !------------------------------------------------------------------------------- IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN iml = assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1), & & SIZE(lon_in2)],TRIM(modname)//" iml") jml = assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2), & & TRIM(modname)//" jml") lml = assert_eq( SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3), & & TRIM(modname)//" lml") jml2 = SIZE(lat_in2) !--- CHECK IF THE FIELD IS KNOWN SELECT CASE(var) CASE('u'); vname='U' CASE('v'); vname='V' CASE('t'); vname='TEMP' CASE('q'); vname='R'; msg='humidity as the saturated humidity' CASE('tpot'); vname='TEMP'; msg='potential temperature as the Exner function' CASE DEFAULT; msg='No rule to extract variable '//TRIM(var) CALL abort_gcm(modname,TRIM(msg)//' from any data set',1) END SELECT !--- CHECK IF SOMETHING IS MISSING IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN msg='Could not compute '//TRIM(msg)//' is missing or constant.' CALL abort_gcm(modname,TRIM(msg),1) END IF !--- INTERPOLATE 3D FIELD IF NEEDED IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, & lat_in2,pls,champ,ibar) !--- COMPUTE THE REQUIRED FILED SELECT CASE(var) CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO champ(iml,:,:)=champ(1,:,:) !--- Eastward wind CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO champ(iml,:,:)=champ(1,:,:) !--- Northward wind CASE('tpot','q') IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Temperature ELSE; champ=champ*.01*workvar !--- Relat. humidity WHERE(champ<0.) champ=1.0E-10 END IF DO il=1,lml xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols champ(:,1 ,il) = xppn champ(:,jml,il) = xpps END DO END SELECT END IF END SUBROUTINE startget_dyn3d ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque) ! !------------------------------------------------------------------------------- USE conf_dat_m, ONLY: conf_dat2d IMPLICIT NONE !=============================================================================== ! Purpose: Compute "phis" just like it would be in start_init_orog. !=============================================================================== ! Arguments: REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml) !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: modname="start_init_orog0" CHARACTER(LEN=256) :: title="RELIEF" INTEGER :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1) REAL :: lev(1), date, dt REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:) REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:,:) !------------------------------------------------------------------------------- iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml") jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml") IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1) IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1) pi=2.0*ASIN(1.0); deg2rad=pi/180.0 IF(ANY(phis/=-99999.)) RETURN !--- phis ALREADY KNOWN !--- HIGH RESOLUTION OROGRAPHY CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid) ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel)) CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,& lev, ttm_tmp, itau, date, dt, fid) ALLOCATE(relief_hi(iml_rel,jml_rel)) CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi) CALL flinclo(fid) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel)) lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel)) CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.) DEALLOCATE(lon_ini,lat_ini) !--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0 WRITE(lunout,*) WRITE(lunout,*)'*** Compute surface geopotential ***' !--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque) phis = phis * 9.81 phis(iml,:) = phis(1,:) DEALLOCATE(relief_hi,lon_rad,lat_rad) END SUBROUTINE start_init_orog0 ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask) ! !=============================================================================== ! Purpose: Extracted from grid_noro to provide geopotential height for dynamics ! without any call to physics subroutines. !=============================================================================== IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp) REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp) REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) REAL, INTENT(INOUT):: mask(:,:) !--- MASK (imar+1,jmar) !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: modname="grid_noro0" REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2) REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2) REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar) REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:) ! dim (imar+1,jmar) REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax) REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax) REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax) LOGICAL :: masque_lu INTEGER :: i, ii, imdp, imar, iext INTEGER :: j, jj, jmdp, jmar, nn REAL :: xpi, zlenx, weighx, xincr, zbordnor, zmeanor, zweinor, zbordest REAL :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue !------------------------------------------------------------------------------- imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp") jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp") imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1 jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar") IF(imar/=iim) CALL abort_gcm(TRIM(modname),'imar/=iim' ,1) IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1) iext=imdp/10 xpi = ACOS(-1.) rad = 6371229. !--- ARE WE USING A READ MASK ? masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 WRITE(lunout,*)'Masque lu: ',masque_lu !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES: ALLOCATE(xusn(imdp+2*iext)) xusn(1 +iext:imdp +iext)=xd(:) xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi ALLOCATE(yusn(jmdp+2)) yusn(1 )=yd(1) +(yd(1) -yd(2)) yusn(2:jmdp+1)=yd(:) yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1)) ALLOCATE(zusn(imdp+2*iext,jmdp+2)) zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :) zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :) zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :) zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2) zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2) zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1) zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1) !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) ALLOCATE(a(imar+1),b(imar+1)) b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0 b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0 a(1)=x(1)-(x(2)-x(1))/2.0 a(2:imar+1)= b(1:imar) ALLOCATE(c(jmar),d(jmar)) d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0 d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0 c(1)=y(1)-(y(2)-y(1))/2.0 c(2:jmar)=d(1:jmar-1) !--- INITIALIZATIONS: ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0 ALLOCATE(zmea (imar+1,jmar)); zmea (:,:)= 0.0 !--- SUMMATION OVER GRIDPOINT AREA zleny=xpi/REAL(jmdp)*rad xincr=xpi/REAL(jmdp)/2. ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0. ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0. DO ii = 1, imar+1 DO jj = 1, jmar DO j = 2,jmdp+1 zlenx =zleny *COS(yusn(j)) zbordnor=(xincr+c(jj)-yusn(j))*rad zbordsud=(xincr-d(jj)+yusn(j))*rad weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny)) IF(weighy/=0) THEN DO i = 2, imdp+2*iext-1 zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j)) zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j)) weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx)) IF(weighx/=0)THEN num_tot(ii,jj)=num_tot(ii,jj)+1.0 IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0 weight(ii,jj)=weight(ii,jj)+weighx*weighy zmea (ii,jj)=zmea (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN END IF END DO END IF END DO END DO END DO !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME IF(.NOT.masque_lu) THEN WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:) END IF nn=COUNT(weight(:,1:jmar-1)==0.0) IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:) !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS) ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0 WHERE(mask>=0.1) mask_tmp = 1. WHERE(weight(:,:)/=0.0) zphi(:,:)=mask_tmp(:,:)*zmea(:,:) zmea(:,:)=mask_tmp(:,:)*zmea(:,:) END WHERE WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea) !--- Values at poles zphi(imar+1,:)=zphi(1,:) zweinor=SUM(weight(1:imar, 1),DIM=1) zweisud=SUM(weight(1:imar,jmar),DIM=1) zmeanor=SUM(weight(1:imar, 1)*zmea(1:imar, 1),DIM=1) zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1) zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud END SUBROUTINE grid_noro0 ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol) ! !------------------------------------------------------------------------------- IMPLICIT NONE !=============================================================================== ! Purpose: Compute psol, knowing phis. !=============================================================================== ! Arguments: REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) LOGICAL, INTENT(IN) :: ibar REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml) REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml) !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: modname='start_init_dyn' REAL :: date, dt INTEGER :: iml, jml, jml2, itau(1) REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:) REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) REAL, ALLOCATABLE :: z(:,:), ps(:,:), ts(:,:) !------------------------------------------------------------------------------- iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2), & & TRIM(modname)//" iml") jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml") jml2=SIZE(lat_in2) WRITE(lunout,*) 'Opening the surface analysis' CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn)) ALLOCATE(levdyn_ini(llm_dyn)) CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, & lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn) !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) CALL get_var_dyn('Z',z) !--- SURFACE GEOPOTENTIAL CALL get_var_dyn('SP',ps) !--- SURFACE PRESSURE CALL get_var_dyn('ST',ts) !--- SURFACE TEMPERATURE ! CALL flinclo(fid_dyn) DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) !--- PSOL IS COMPUTED IN PASCALS psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0 & & /ts(:iml-1,:)) psol(iml,:)=psol(1,:) DEALLOCATE(z,ps,ts) psol(:,1 )=SUM(aire(1:iml-1,1 )*psol(1:iml-1,1 ))/apoln !--- NORTH POLE psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols !--- SOUTH POLE CONTAINS !------------------------------------------------------------------------------- ! SUBROUTINE get_var_dyn(title,field) ! !------------------------------------------------------------------------------- USE conf_dat_m, ONLY: conf_dat2d IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: title REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:) !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: msg INTEGER :: tllm !------------------------------------------------------------------------------- SELECT CASE(title) CASE('Z'); tllm=0; msg='geopotential' CASE('SP'); tllm=0; msg='surface pressure' CASE('ST'); tllm=llm_dyn; msg='temperature' END SELECT IF(.NOT.ALLOCATED(field)) THEN ALLOCATE(field(iml,jml)) CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar) CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana, & lon_in, lat_in, lon_in2, lat_in2, field) ELSE IF(SIZE(field)/=SIZE(z)) THEN msg='The '//TRIM(msg)//' field we have does not have the right size' CALL abort_gcm(TRIM(modname),msg,1) END IF END SUBROUTINE get_var_dyn ! !------------------------------------------------------------------------------- END SUBROUTINE start_init_dyn ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar) ! !------------------------------------------------------------------------------- USE conf_dat_m, ONLY: conf_dat3d USE pchsp_95_m, ONLY: pchsp_95 USE pchfe_95_m, ONLY: pchfe_95 IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: var REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml) REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml) LOGICAL, INTENT(IN) :: ibar !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: modname='start_inter_3d' LOGICAL :: skip REAL :: chmin, chmax INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr INTEGER :: n_extrap ! Extrapolated points number REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:) REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:) REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:) !------------------------------------------------------------------------------- iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml") jml=assert_eq(SIZE(lat_in), SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml") lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2) WRITE(lunout, *)'Going into flinget to extract the 3D field.' IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) !--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & lon_rad, lat_rad, lev_dyn, var_ana3d, ibar) DEALLOCATE(lon_ini, lat_ini) !--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) DO il = 1,llm_dyn CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), & lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il)) END DO DEALLOCATE(lon_rad, lat_rad) !--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) ax = lev_dyn(llm_dyn:1:-1) skip = .FALSE. n_extrap = 0 DO ij=1, jml DO ii=1, iml-1 ay = var_tmp3d(ii, ij, llm_dyn:1:-1) yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.) CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & var3d(ii, ij, lml:1:-1), ierr) IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1) n_extrap = n_extrap + ierr END DO END DO IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap var3d(iml, :, :) = var3d(1, :, :) DO il=1, lml CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax) WRITE(lunout, *)' '//TRIM(var)//' min max l ', il, chmin, chmax END DO END SUBROUTINE start_inter_3d ! !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) ! !------------------------------------------------------------------------------- USE inter_barxy_m, ONLY: inter_barxy USE grid_atob_m, ONLY: grille_m IMPLICIT NONE !------------------------------------------------------------------------------- ! Arguments: CHARACTER(LEN=*), INTENT(IN) :: nam LOGICAL, INTENT(IN) :: ibar, ibeg REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) !------------------------------------------------------------------------------- ! Local variables: CHARACTER(LEN=256) :: modname="interp_startvar" INTEGER :: ii, jj, i1, j1, j2 REAL, ALLOCATABLE :: vtmp(:,:) !------------------------------------------------------------------------------- ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1") j2=SIZE(lat2) ALLOCATE(vtmp(i1-1,j1)) IF(ibar) THEN IF(ibeg.AND.prt_level>1) THEN WRITE(lunout,*)"---------------------------------------------------------" WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" WRITE(lunout,*)"---------------------------------------------------------" END IF CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) ELSE CALL grille_m (lon, lat, vari, lon1, lat1, vtmp) END IF CALL gr_int_dyn(vtmp, varo, i1-1, j1) END SUBROUTINE interp_startvar ! !------------------------------------------------------------------------------- !#endif ! of #ifdef CPP_EARTH END MODULE etat0dyn ! !*******************************************************************************