Changeset 5116 for LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd
- Timestamp:
- Jul 24, 2024, 2:54:37 PM (6 months ago)
- Location:
- LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd
- Files:
-
- 6 edited
-
callphysiq_mod.F90 (modified) (2 diffs)
-
ce0l.F90 (modified) (1 diff)
-
etat0dyn_netcdf.F90 (modified) (2 diffs)
-
etat0phys_netcdf.F90 (modified) (3 diffs)
-
limit_netcdf.f90 (modified) (1 diff)
-
test_disvert_m.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/callphysiq_mod.F90
r5110 r5116 69 69 70 70 !$OMP MASTER 71 if (ok_dyn_xios) then71 if (ok_dyn_xios) THEN 72 72 CALL xios_get_current_context(dyn3d_ctx_handle) 73 73 endif … … 99 99 ! switching back to LMDZDYN context 100 100 !$OMP MASTER 101 if (ok_dyn_xios) then101 if (ok_dyn_xios) THEN 102 102 CALL xios_set_current_context(dyn3d_ctx_handle) 103 103 endif -
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/ce0l.F90
r5106 r5116 36 36 USE mod_hallo, ONLY: init_mod_hallo 37 37 USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys 38 USE lmdz_xios, only: using_xios, xios_finalize38 USE lmdz_xios, ONLY: using_xios, xios_finalize 39 39 #endif 40 40 -
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90
r5113 r5116 1 1 MODULE etat0dyn 2 2 3 !*******************************************************************************4 ! Purpose: Create dynamical initial state using atmospheric fields from a5 ! database of atmospheric to initialize the model.6 !-------------------------------------------------------------------------------7 ! Comments:8 9 ! * This module is designed to work for Earth (and with ioipsl)10 11 ! * etat0dyn_netcdf routine can access to NetCDF data through the following12 ! routine (to be called after restget):13 ! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,&14 ! champ, lon_in2, lat_in2)15 16 ! * Variables should have the following names in the NetCDF files:17 ! 'U' : East ward wind (in "ECDYN.nc")18 ! 'V' : Northward wind (in "ECDYN.nc")19 ! 'TEMP' : Temperature (in "ECDYN.nc")20 ! 'R' : Relative humidity (in "ECDYN.nc")21 ! 'RELIEF' : High resolution orography (in "Relief.nc") 22 23 ! * The land mask and corresponding weights can be:24 ! 1) already known (in particular if etat0dyn has been called before) ;25 ! in this case, ANY(masque(:,:)/=-99999.) = .TRUE.26 ! 2) computed using the ocean mask from the ocean model (to ensure ocean27 ! fractions are the same for atmosphere and ocean) for coupled runs.28 ! File name: "o2a.nc" ; variable name: "OceMask"29 ! 3) computed from topography file "Relief.nc" for forced runs.30 31 ! * There is a big mess with the longitude size. Should it be iml or iml+1 ?32 ! I have chosen to use the iml+1 as an argument to this routine and we declare33 ! internaly smaller fields when needed. This needs to be cleared once and for34 ! all in LMDZ. A convention is required.35 !-------------------------------------------------------------------------------36 USE ioipsl, 37 USE lmdz_assert_eq, 3 !******************************************************************************* 4 ! Purpose: Create dynamical initial state using atmospheric fields from a 5 ! database of atmospheric to initialize the model. 6 !------------------------------------------------------------------------------- 7 ! Comments: 8 9 ! * This module is designed to work for Earth (and with ioipsl) 10 11 ! * etat0dyn_netcdf routine can access to NetCDF data through the following 12 ! routine (to be called after restget): 13 ! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,& 14 ! champ, lon_in2, lat_in2) 15 16 ! * Variables should have the following names in the NetCDF files: 17 ! 'U' : East ward wind (in "ECDYN.nc") 18 ! 'V' : Northward wind (in "ECDYN.nc") 19 ! 'TEMP' : Temperature (in "ECDYN.nc") 20 ! 'R' : Relative humidity (in "ECDYN.nc") 21 ! 'RELIEF' : High resolution orography (in "Relief.nc") 22 23 ! * The land mask and corresponding weights can be: 24 ! 1) already known (in particular if etat0dyn has been called before) ; 25 ! in this case, ANY(masque(:,:)/=-99999.) = .TRUE. 26 ! 2) computed using the ocean mask from the ocean model (to ensure ocean 27 ! fractions are the same for atmosphere and ocean) for coupled runs. 28 ! File name: "o2a.nc" ; variable name: "OceMask" 29 ! 3) computed from topography file "Relief.nc" for forced runs. 30 31 ! * There is a big mess with the longitude size. Should it be iml or iml+1 ? 32 ! I have chosen to use the iml+1 as an argument to this routine and we declare 33 ! internaly smaller fields when needed. This needs to be cleared once and for 34 ! all in LMDZ. A convention is required. 35 !------------------------------------------------------------------------------- 36 USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo 37 USE lmdz_assert_eq, ONLY: assert_eq 38 38 USE comconst_mod, ONLY: pi, cpp, kappa 39 39 USE comvert_mod, ONLY: ap, bp, preff, pressure_exner 40 40 USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itau_phy, start_time 41 41 USE strings_mod, ONLY: strLower 42 42 43 43 IMPLICIT NONE 44 44 … … 52 52 include "comdissnew.h" 53 53 REAL, SAVE :: deg2rad 54 INTEGER, SAVE:: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn55 REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)56 CHARACTER(LEN =120), PARAMETER :: dynfname='ECDYN.nc'54 INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn 55 REAL, ALLOCATABLE, SAVE :: lon_dyn(:, :), lat_dyn(:, :), levdyn_ini(:) 56 CHARACTER(LEN = 120), PARAMETER :: dynfname = 'ECDYN.nc' 57 57 58 58 CONTAINS 59 59 60 !------------------------------------------------------------------------------- 61 62 SUBROUTINE etat0dyn_netcdf(masque, phis) 63 64 !------------------------------------------------------------------------------- 65 ! Purpose: Create dynamical initial states. 66 !------------------------------------------------------------------------------- 67 ! Notes: 1) This routine is designed to work for Earth 68 ! 2) If masque(:,:)/=-99999., masque and phis are already known. 69 ! Otherwise: compute it. 70 !------------------------------------------------------------------------------- 71 USE control_mod 72 USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz 73 USE regr_pr_o3_m, ONLY: regr_pr_o3 74 USE press_coefoz_m, ONLY: press_coefoz 75 USE exner_hyb_m, ONLY: exner_hyb 76 USE exner_milieu_m, ONLY: exner_milieu 77 USE infotrac, ONLY: nqtot, tracers 78 USE lmdz_filtreg 79 USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA 80 IMPLICIT NONE 81 !------------------------------------------------------------------------------- 82 ! Arguments: 83 REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask 84 REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential 85 !------------------------------------------------------------------------------- 86 ! Local variables: 87 CHARACTER(LEN=256) :: modname, fmt 88 INTEGER :: i, j, l, ji, itau, iday, iq 89 REAL :: xpn, xps, time, phystep 90 REAL, DIMENSION(iip1,jjp1) :: psol 91 REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d 92 REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd 93 REAL, DIMENSION(iip1,jjp1,llm) :: pk, pls, y, masse 94 REAL, DIMENSION(iip1,jjm ,llm) :: vvent 95 REAL, DIMENSION(ip1jm ,llm) :: pbarv 96 REAL, DIMENSION(ip1jmp1 ,llm) :: pbaru, phi, w 97 REAL, DIMENSION(ip1jmp1) :: pks 98 REAL, DIMENSION(iim) :: xppn, xpps 99 REAL, ALLOCATABLE :: q3d(:,:,:,:) 100 !------------------------------------------------------------------------------- 101 modname='etat0dyn_netcdf' 102 103 deg2rad = pi/180.0 104 y(:,:,:)=0 !ym warning unitialized variable 105 106 ! Compute psol AND tsol, knowing phis. 107 !******************************************************************************* 108 CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol) 109 110 ! Mid-levels pressure computation 111 !******************************************************************************* 112 CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d 113 IF(pressure_exner) THEN !--- Update pk, pks 114 CALL exner_hyb (ip1jmp1,psol,p3d,pks,pk) 115 ELSE 116 CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk) 117 END IF 118 pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa) !--- Update pls 119 120 ! Update uvent, vvent, t3d and tpot 121 !******************************************************************************* 122 uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0 123 CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv) 124 CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, & 125 rlonu,rlatu(:jjm)) 126 CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv) 127 tpot(:,:,:)=t3d(:,:,:) 128 CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv) 129 130 WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:)) 131 WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:)) 132 133 ! Humidity at saturation computation 134 !******************************************************************************* 135 WRITE(lunout,*) 'avant q_sat' 136 CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat) 137 WRITE(lunout,*) 'apres q_sat' 138 WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:)) 139 ! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) 140 qd (:,:,:) = 0.0 141 CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv) 142 ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:) 143 CALL flinclo(fid_dyn) 144 145 ! Parameterization of ozone chemistry: 146 !******************************************************************************* 147 ! Look for ozone tracer: 148 IF (CPPKEY_INCA) THEN 149 DO iq=1,nqtot; IF(strLower(tracers(iq)%name)=="o3") EXIT; END DO 150 IF(iq/=nqtot+1) THEN 151 CALL regr_lat_time_coefoz 152 CALL press_coefoz 153 CALL regr_pr_o3(p3d, q3d(:,:,:,iq)) 154 q3d(:,:,:,iq)=q3d(:,:,:,iq)*48./ 29. !--- Mole->mass fraction 155 END IF 156 END IF 157 q3d(iip1,:,:,:)=q3d(1,:,:,:) 158 159 ! Writing 160 !******************************************************************************* 161 CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & 162 tetatemp, vert_prof_dissip) 163 WRITE(lunout,*)'sortie inidissip' 164 itau=0 165 itau_dyn=0 166 itau_phy=0 167 iday=dayref+itau/day_step 168 time=FLOAT(itau-(iday-dayref)*day_step)/day_step 169 IF(time>1.) THEN 170 time=time-1 171 iday=iday+1 172 END IF 173 day_ref=dayref 174 annee_ref=anneeref 175 CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) 176 WRITE(lunout,*)'sortie geopot' 177 CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & 178 phi, w, pbaru, pbarv, time+iday-dayref) 179 WRITE(lunout,*)'sortie caldyn0' 180 start_time = 0. 60 !------------------------------------------------------------------------------- 61 62 SUBROUTINE etat0dyn_netcdf(masque, phis) 63 64 !------------------------------------------------------------------------------- 65 ! Purpose: Create dynamical initial states. 66 !------------------------------------------------------------------------------- 67 ! Notes: 1) This routine is designed to work for Earth 68 ! 2) If masque(:,:)/=-99999., masque and phis are already known. 69 ! Otherwise: compute it. 70 !------------------------------------------------------------------------------- 71 USE control_mod 72 USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz 73 USE regr_pr_o3_m, ONLY: regr_pr_o3 74 USE press_coefoz_m, ONLY: press_coefoz 75 USE exner_hyb_m, ONLY: exner_hyb 76 USE exner_milieu_m, ONLY: exner_milieu 77 USE infotrac, ONLY: nqtot, tracers 78 USE lmdz_filtreg 79 USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA 80 IMPLICIT NONE 81 !------------------------------------------------------------------------------- 82 ! Arguments: 83 REAL, INTENT(INOUT) :: masque(iip1, jjp1) !--- Land-ocean mask 84 REAL, INTENT(INOUT) :: phis (iip1, jjp1) !--- Ground geopotential 85 !------------------------------------------------------------------------------- 86 ! Local variables: 87 CHARACTER(LEN = 256) :: modname, fmt 88 INTEGER :: i, j, l, ji, itau, iday, iq 89 REAL :: xpn, xps, time, phystep 90 REAL, DIMENSION(iip1, jjp1) :: psol 91 REAL, DIMENSION(iip1, jjp1, llm + 1) :: p3d 92 REAL, DIMENSION(iip1, jjp1, llm) :: uvent, t3d, tpot, qsat, qd 93 REAL, DIMENSION(iip1, jjp1, llm) :: pk, pls, y, masse 94 REAL, DIMENSION(iip1, jjm, llm) :: vvent 95 REAL, DIMENSION(ip1jm, llm) :: pbarv 96 REAL, DIMENSION(ip1jmp1, llm) :: pbaru, phi, w 97 REAL, DIMENSION(ip1jmp1) :: pks 98 REAL, DIMENSION(iim) :: xppn, xpps 99 REAL, ALLOCATABLE :: q3d(:, :, :, :) 100 !------------------------------------------------------------------------------- 101 modname = 'etat0dyn_netcdf' 102 103 deg2rad = pi / 180.0 104 y(:, :, :) = 0 !ym warning unitialized variable 105 106 ! Compute psol AND tsol, knowing phis. 107 !******************************************************************************* 108 CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol) 109 110 ! Mid-levels pressure computation 111 !******************************************************************************* 112 CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d 113 IF(pressure_exner) THEN !--- Update pk, pks 114 CALL exner_hyb (ip1jmp1, psol, p3d, pks, pk) 115 ELSE 116 CALL exner_milieu(ip1jmp1, psol, p3d, pks, pk) 117 END IF 118 pls(:, :, :) = preff * (pk(:, :, :) / cpp)**(1. / kappa) !--- Update pls 119 120 ! Update uvent, vvent, t3d and tpot 121 !******************************************************************************* 122 uvent(:, :, :) = 0.0 ; vvent(:, :, :) = 0.0 ; t3d (:, :, :) = 0.0 123 CALL startget_dyn3d('u', rlonu, rlatu, pls, y, uvent, rlonv, rlatv) 124 CALL startget_dyn3d('v', rlonv, rlatv, pls(:, :jjm, :), y(:, :jjm, :), vvent, & 125 rlonu, rlatu(:jjm)) 126 CALL startget_dyn3d('t', rlonv, rlatu, pls, y, t3d, rlonu, rlatv) 127 tpot(:, :, :) = t3d(:, :, :) 128 CALL startget_dyn3d('tpot', rlonv, rlatu, pls, pk, tpot, rlonu, rlatv) 129 130 WRITE(lunout, *) 'T3D min,max:', MINVAL(t3d(:, :, :)), MAXVAL(t3d(:, :, :)) 131 WRITE(lunout, *) 'PLS min,max:', MINVAL(pls(:, :, :)), MAXVAL(pls(:, :, :)) 132 133 ! Humidity at saturation computation 134 !******************************************************************************* 135 WRITE(lunout, *) 'avant q_sat' 136 CALL q_sat(llm * jjp1 * iip1, t3d, pls, qsat) 137 WRITE(lunout, *) 'apres q_sat' 138 WRITE(lunout, *) 'QSAT min,max:', MINVAL(qsat(:, :, :)), MAXVAL(qsat(:, :, :)) 139 ! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) 140 qd (:, :, :) = 0.0 141 CALL startget_dyn3d('q', rlonv, rlatu, pls, qsat, qd, rlonu, rlatv) 142 ALLOCATE(q3d(iip1, jjp1, llm, nqtot)); q3d(:, :, :, :) = 0.0 ; q3d(:, :, :, 1) = qd(:, :, :) 143 CALL flinclo(fid_dyn) 144 145 ! Parameterization of ozone chemistry: 146 !******************************************************************************* 147 ! Look for ozone tracer: 148 IF (CPPKEY_INCA) THEN 149 DO iq = 1, nqtot; IF(strLower(tracers(iq)%name)=="o3") EXIT; 150 END DO 151 IF(iq/=nqtot + 1) THEN 152 CALL regr_lat_time_coefoz 153 CALL press_coefoz 154 CALL regr_pr_o3(p3d, q3d(:, :, :, iq)) 155 q3d(:, :, :, iq) = q3d(:, :, :, iq) * 48. / 29. !--- Mole->mass fraction 156 END IF 157 END IF 158 q3d(iip1, :, :, :) = q3d(1, :, :, :) 159 160 ! Writing 161 !******************************************************************************* 162 CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & 163 tetatemp, vert_prof_dissip) 164 WRITE(lunout, *)'sortie inidissip' 165 itau = 0 166 itau_dyn = 0 167 itau_phy = 0 168 iday = dayref + itau / day_step 169 time = FLOAT(itau - (iday - dayref) * day_step) / day_step 170 IF(time>1.) THEN 171 time = time - 1 172 iday = iday + 1 173 END IF 174 day_ref = dayref 175 annee_ref = anneeref 176 CALL geopot(ip1jmp1, tpot, pk, pks, phis, phi) 177 WRITE(lunout, *)'sortie geopot' 178 CALL caldyn0(itau, uvent, vvent, tpot, psol, masse, pk, phis, & 179 phi, w, pbaru, pbarv, time + iday - dayref) 180 WRITE(lunout, *)'sortie caldyn0' 181 start_time = 0. 181 182 #ifdef CPP_PARA 182 183 CALL dynredem0_loc( "start.nc", dayref, phis) 183 184 #else 184 CALL dynredem0("start.nc", dayref, phis)185 CALL dynredem0("start.nc", dayref, phis) 185 186 #endif 186 WRITE(lunout,*)'sortie dynredem0'187 WRITE(lunout, *)'sortie dynredem0' 187 188 #ifdef CPP_PARA 188 189 CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 189 190 #else 190 CALL dynredem1("start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)191 CALL dynredem1("start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) 191 192 #endif 192 WRITE(lunout,*)'sortie dynredem1' 193 CALL histclo() 194 195 END SUBROUTINE etat0dyn_netcdf 196 197 !------------------------------------------------------------------------------- 198 199 200 !------------------------------------------------------------------------------- 201 202 SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,& 203 champ, lon_in2, lat_in2) 204 !------------------------------------------------------------------------------- 205 IMPLICIT NONE 206 !=============================================================================== 207 ! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R 208 ! (3D fields) of file dynfname. 209 !------------------------------------------------------------------------------- 210 ! Note: An input auxilliary field "workvar" has to be specified in two cases: 211 ! * for "q": the saturated humidity. 212 ! * for "tpot": the Exner function. 213 !=============================================================================== 214 ! Arguments: 215 CHARACTER(LEN=*), INTENT(IN) :: var 216 REAL, INTENT(IN) :: lon_in(:) ! dim (iml) 217 REAL, INTENT(IN) :: lat_in(:) ! dim (jml) 218 REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml) 219 REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml) 220 REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml) 221 REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) 222 REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) 223 !------------------------------------------------------------------------------- 224 ! Local variables: 225 CHARACTER(LEN=10) :: vname 226 CHARACTER(LEN=256) :: msg, modname="startget_dyn3d" 227 INTEGER :: iml, jml, jml2, lml, il 228 REAL :: xppn, xpps 229 !------------------------------------------------------------------------------- 230 iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1), & 231 SIZE(lon_in2)], TRIM(modname)//" iml") 232 jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2), & 233 TRIM(modname)//" jml") 234 lml=assert_eq( SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3), & 235 TRIM(modname)//" lml") 236 jml2=SIZE(lat_in2) 237 238 !--- CHECK IF THE FIELD IS KNOWN 239 SELECT CASE(var) 240 CASE('u'); vname='U' 241 CASE('v'); vname='V' 242 CASE('t'); vname='TEMP' 243 CASE('q'); vname='R'; msg='humidity as the saturated humidity' 244 CASE('tpot'); msg='potential temperature as the Exner function' 245 CASE DEFAULT; msg='No rule to extract variable '//TRIM(var) 246 CALL abort_gcm(modname,TRIM(msg)//' from any data set',1) 247 END SELECT 248 249 !--- CHECK IF SOMETHING IS MISSING 250 IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN 251 msg='Could not compute '//TRIM(msg)//' is missing or constant.' 252 CALL abort_gcm(modname,TRIM(msg),1) 253 END IF 254 255 !--- INTERPOLATE 3D FIELD IF NEEDED 256 IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, & 257 lat_in2,pls,champ) 258 259 !--- COMPUTE THE REQUIRED FILED 260 SELECT CASE(var) 261 CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO 262 champ(iml,:,:)=champ(1,:,:) !--- Eastward wind 263 264 CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO 265 champ(iml,:,:)=champ(1,:,:) !--- Northward wind 266 267 CASE('tpot','q') 268 IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature 269 ELSE; champ=champ*.01*workvar !--- Relative humidity 270 WHERE(champ<0.) champ=1.0E-10 193 WRITE(lunout, *)'sortie dynredem1' 194 CALL histclo() 195 196 END SUBROUTINE etat0dyn_netcdf 197 198 !------------------------------------------------------------------------------- 199 200 201 !------------------------------------------------------------------------------- 202 203 SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar, & 204 champ, lon_in2, lat_in2) 205 !------------------------------------------------------------------------------- 206 IMPLICIT NONE 207 !=============================================================================== 208 ! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R 209 ! (3D fields) of file dynfname. 210 !------------------------------------------------------------------------------- 211 ! Note: An input auxilliary field "workvar" has to be specified in two cases: 212 ! * for "q": the saturated humidity. 213 ! * for "tpot": the Exner function. 214 !=============================================================================== 215 ! Arguments: 216 CHARACTER(LEN = *), INTENT(IN) :: var 217 REAL, INTENT(IN) :: lon_in(:) ! dim (iml) 218 REAL, INTENT(IN) :: lat_in(:) ! dim (jml) 219 REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml) 220 REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml) 221 REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml) 222 REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) 223 REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) 224 !------------------------------------------------------------------------------- 225 ! Local variables: 226 CHARACTER(LEN = 10) :: vname 227 CHARACTER(LEN = 256) :: msg, modname = "startget_dyn3d" 228 INTEGER :: iml, jml, jml2, lml, il 229 REAL :: xppn, xpps 230 !------------------------------------------------------------------------------- 231 iml = assert_eq([SIZE(lon_in), SIZE(pls, 1), SIZE(workvar, 1), SIZE(champ, 1), & 232 SIZE(lon_in2)], TRIM(modname) // " iml") 233 jml = assert_eq(SIZE(lat_in), SIZE(pls, 2), SIZE(workvar, 2), SIZE(champ, 2), & 234 TRIM(modname) // " jml") 235 lml = assert_eq(SIZE(pls, 3), SIZE(workvar, 3), SIZE(champ, 3), & 236 TRIM(modname) // " lml") 237 jml2 = SIZE(lat_in2) 238 239 !--- CHECK IF THE FIELD IS KNOWN 240 SELECT CASE(var) 241 CASE('u'); vname = 'U' 242 CASE('v'); vname = 'V' 243 CASE('t'); vname = 'TEMP' 244 CASE('q'); vname = 'R'; msg = 'humidity as the saturated humidity' 245 CASE('tpot'); msg = 'potential temperature as the Exner function' 246 CASE DEFAULT; msg = 'No rule to extract variable ' // TRIM(var) 247 CALL abort_gcm(modname, TRIM(msg) // ' from any data set', 1) 248 END SELECT 249 250 !--- CHECK IF SOMETHING IS MISSING 251 IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN 252 msg = 'Could not compute ' // TRIM(msg) // ' is missing or constant.' 253 CALL abort_gcm(modname, TRIM(msg), 1) 254 END IF 255 256 !--- INTERPOLATE 3D FIELD IF NEEDED 257 IF(var/='tpot') CALL start_inter_3d(TRIM(vname), lon_in, lat_in, lon_in2, & 258 lat_in2, pls, champ) 259 260 !--- COMPUTE THE REQUIRED FILED 261 SELECT CASE(var) 262 CASE('u'); DO il = 1, lml; champ(:, :, il) = champ(:, :, il) * cu(:, 1:jml); 263 END DO 264 champ(iml, :, :) = champ(1, :, :) !--- Eastward wind 265 266 CASE('v'); DO il = 1, lml; champ(:, :, il) = champ(:, :, il) * cv(:, 1:jml); 267 END DO 268 champ(iml, :, :) = champ(1, :, :) !--- Northward wind 269 270 CASE('tpot', 'q') 271 IF(var=='tpot') THEN; champ = champ * cpp / workvar !--- Potential temperature 272 ELSE; champ = champ * .01 * workvar !--- Relative humidity 273 WHERE(champ<0.) champ = 1.0E-10 271 274 END IF 272 DO il =1,lml273 xppn = SUM(aire(:, 1 )*champ(:,1 ,il))/apoln274 xpps = SUM(aire(:, jml)*champ(:,jml,il))/apols275 champ(:, 1 ,il) = xppn276 champ(:, jml,il) = xpps275 DO il = 1, lml 276 xppn = SUM(aire(:, 1) * champ(:, 1, il)) / apoln 277 xpps = SUM(aire(:, jml) * champ(:, jml, il)) / apols 278 champ(:, 1, il) = xppn 279 champ(:, jml, il) = xpps 277 280 END DO 278 END SELECT 279 280 END SUBROUTINE startget_dyn3d 281 282 !------------------------------------------------------------------------------- 283 284 285 !------------------------------------------------------------------------------- 286 287 SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol) 288 289 !------------------------------------------------------------------------------- 290 IMPLICIT NONE 291 !=============================================================================== 292 ! Purpose: Compute psol, knowing phis. 293 !=============================================================================== 294 ! Arguments: 295 REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) 296 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 297 REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml) 298 REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml) 299 !------------------------------------------------------------------------------- 300 ! Local variables: 301 CHARACTER(LEN=256) :: modname='start_init_dyn' 302 REAL :: date, dt 303 INTEGER :: iml, jml, jml2, itau(1) 304 REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:) 305 REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) 306 REAL, ALLOCATABLE :: z(:,:), ps(:,:), ts(:,:) 307 !------------------------------------------------------------------------------- 308 iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2), & 309 TRIM(modname)//" iml") 310 jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml") 311 jml2=SIZE(lat_in2) 312 313 WRITE(lunout,*) 'Opening the surface analysis' 314 CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) 315 WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn 316 317 ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn)) 318 ALLOCATE(levdyn_ini(llm_dyn)) 319 CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, & 320 lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn) 321 322 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 323 ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) 324 lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad 325 lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad 326 327 ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) 328 CALL get_var_dyn('Z',z) !--- SURFACE GEOPOTENTIAL 329 CALL get_var_dyn('SP',ps) !--- SURFACE PRESSURE 330 CALL get_var_dyn('ST',ts) !--- SURFACE TEMPERATURE 331 ! CALL flinclo(fid_dyn) 332 DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) 333 334 !--- PSOL IS COMPUTED IN PASCALS 335 psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0 & 336 /ts(:iml-1,:)) 337 psol(iml,:)=psol(1,:) 338 DEALLOCATE(z,ps,ts) 339 psol(:,1 )=SUM(aire(1:iml-1,1 )*psol(1:iml-1,1 ))/apoln !--- NORTH POLE 340 psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols !--- SOUTH POLE 341 342 CONTAINS 343 344 !------------------------------------------------------------------------------- 345 346 SUBROUTINE get_var_dyn(title,field) 347 348 !------------------------------------------------------------------------------- 349 USE conf_dat_m, ONLY: conf_dat2d 350 IMPLICIT NONE 351 !------------------------------------------------------------------------------- 352 ! Arguments: 353 CHARACTER(LEN=*), INTENT(IN) :: title 354 REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:) 355 !------------------------------------------------------------------------------- 356 ! Local variables: 357 CHARACTER(LEN=256) :: msg 358 INTEGER :: tllm 359 !------------------------------------------------------------------------------- 360 SELECT CASE(title) 361 CASE('Z'); tllm=0; msg='geopotential' 362 CASE('SP'); tllm=0; msg='surface pressure' 363 CASE('ST'); tllm=llm_dyn; msg='temperature' 364 END SELECT 365 IF(.NOT.ALLOCATED(field)) THEN 366 ALLOCATE(field(iml,jml)) 367 CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) 368 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) 369 CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana, & 370 lon_in, lat_in, lon_in2, lat_in2, field) 371 ELSE IF(SIZE(field)/=SIZE(z)) THEN 372 msg='The '//TRIM(msg)//' field we have does not have the right size' 373 CALL abort_gcm(TRIM(modname),msg,1) 374 END IF 375 376 END SUBROUTINE get_var_dyn 377 378 !------------------------------------------------------------------------------- 379 380 END SUBROUTINE start_init_dyn 381 382 !------------------------------------------------------------------------------- 383 384 385 !------------------------------------------------------------------------------- 386 387 SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d) 388 389 !------------------------------------------------------------------------------- 390 USE conf_dat_m, ONLY: conf_dat3d 391 USE lmdz_libmath_pch, ONLY: pchsp_95, pchfe_95 392 IMPLICIT NONE 393 !------------------------------------------------------------------------------- 394 ! Arguments: 395 CHARACTER(LEN=*), INTENT(IN) :: var 396 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) 397 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 398 REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml) 399 REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml) 400 !------------------------------------------------------------------------------- 401 ! Local variables: 402 CHARACTER(LEN=256) :: modname='start_inter_3d' 403 LOGICAL :: skip 404 REAL :: chmin, chmax 405 INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr 406 INTEGER :: n_extrap ! Extrapolated points number 407 REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:) 408 REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:) 409 REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:) 410 !------------------------------------------------------------------------------- 411 iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml") 412 jml=assert_eq(SIZE(lat_in), SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml") 413 lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2) 414 415 WRITE(lunout, *)'Going into flinget to extract the 3D field.' 416 IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) 417 CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) 418 419 !--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS 420 ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) 421 lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad 422 lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad 423 424 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 425 ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) 426 CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & 427 lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.) 428 DEALLOCATE(lon_ini, lat_ini) 429 430 !--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro 431 ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) 432 DO il = 1,llm_dyn 433 CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), & 434 lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il)) 435 END DO 436 DEALLOCATE(lon_rad, lat_rad) 437 438 !--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND 439 ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) 440 ax = lev_dyn(llm_dyn:1:-1) 441 skip = .FALSE. 442 n_extrap = 0 443 DO ij=1, jml 444 DO ii=1, iml-1 445 ay = var_tmp3d(ii, ij, llm_dyn:1:-1) 446 yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.) 447 CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & 448 var3d(ii, ij, lml:1:-1), ierr) 449 IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1) 450 n_extrap = n_extrap + ierr 281 END SELECT 282 283 END SUBROUTINE startget_dyn3d 284 285 !------------------------------------------------------------------------------- 286 287 288 !------------------------------------------------------------------------------- 289 290 SUBROUTINE start_init_dyn(lon_in, lat_in, lon_in2, lat_in2, zs, psol) 291 292 !------------------------------------------------------------------------------- 293 IMPLICIT NONE 294 !=============================================================================== 295 ! Purpose: Compute psol, knowing phis. 296 !=============================================================================== 297 ! Arguments: 298 REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) 299 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 300 REAL, INTENT(IN) :: zs (:, :) ! dim (iml,jml) 301 REAL, INTENT(OUT) :: psol(:, :) ! dim (iml,jml) 302 !------------------------------------------------------------------------------- 303 ! Local variables: 304 CHARACTER(LEN = 256) :: modname = 'start_init_dyn' 305 REAL :: date, dt 306 INTEGER :: iml, jml, jml2, itau(1) 307 REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:, :) 308 REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) 309 REAL, ALLOCATABLE :: z(:, :), ps(:, :), ts(:, :) 310 !------------------------------------------------------------------------------- 311 iml = assert_eq(SIZE(lon_in), SIZE(zs, 1), SIZE(psol, 1), SIZE(lon_in2), & 312 TRIM(modname) // " iml") 313 jml = assert_eq(SIZE(lat_in), SIZE(zs, 2), SIZE(psol, 2), TRIM(modname) // " jml") 314 jml2 = SIZE(lat_in2) 315 316 WRITE(lunout, *) 'Opening the surface analysis' 317 CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) 318 WRITE(lunout, *) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn 319 320 ALLOCATE(lon_dyn(iml_dyn, jml_dyn), lat_dyn(iml_dyn, jml_dyn)) 321 ALLOCATE(levdyn_ini(llm_dyn)) 322 CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, & 323 lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn) 324 325 !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS 326 ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) 327 lon_ini(:) = lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini = lon_ini * deg2rad 328 lat_ini(:) = lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini = lat_ini * deg2rad 329 330 ALLOCATE(var_ana(iml_dyn, jml_dyn), lon_rad(iml_dyn), lat_rad(jml_dyn)) 331 CALL get_var_dyn('Z', z) !--- SURFACE GEOPOTENTIAL 332 CALL get_var_dyn('SP', ps) !--- SURFACE PRESSURE 333 CALL get_var_dyn('ST', ts) !--- SURFACE TEMPERATURE 334 ! CALL flinclo(fid_dyn) 335 DEALLOCATE(var_ana, lon_rad, lat_rad, lon_ini, lat_ini) 336 337 !--- PSOL IS COMPUTED IN PASCALS 338 psol(:iml - 1, :) = ps(:iml - 1, :) * (1.0 + (z(:iml - 1, :) - zs(:iml - 1, :)) / 287.0 & 339 / ts(:iml - 1, :)) 340 psol(iml, :) = psol(1, :) 341 DEALLOCATE(z, ps, ts) 342 psol(:, 1) = SUM(aire(1:iml - 1, 1) * psol(1:iml - 1, 1)) / apoln !--- NORTH POLE 343 psol(:, jml) = SUM(aire(1:iml - 1, jml) * psol(1:iml - 1, jml)) / apols !--- SOUTH POLE 344 345 CONTAINS 346 347 !------------------------------------------------------------------------------- 348 349 SUBROUTINE get_var_dyn(title, field) 350 351 !------------------------------------------------------------------------------- 352 USE conf_dat_m, ONLY: conf_dat2d 353 IMPLICIT NONE 354 !------------------------------------------------------------------------------- 355 ! Arguments: 356 CHARACTER(LEN = *), INTENT(IN) :: title 357 REAL, ALLOCATABLE, INTENT(INOUT) :: field(:, :) 358 !------------------------------------------------------------------------------- 359 ! Local variables: 360 CHARACTER(LEN = 256) :: msg 361 INTEGER :: tllm 362 !------------------------------------------------------------------------------- 363 SELECT CASE(title) 364 CASE('Z'); tllm = 0; msg = 'geopotential' 365 CASE('SP'); tllm = 0; msg = 'surface pressure' 366 CASE('ST'); tllm = llm_dyn; msg = 'temperature' 367 END SELECT 368 IF(.NOT.ALLOCATED(field)) THEN 369 ALLOCATE(field(iml, jml)) 370 CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) 371 CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) 372 CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, & 373 lon_in, lat_in, lon_in2, lat_in2, field) 374 ELSE IF(SIZE(field)/=SIZE(z)) THEN 375 msg = 'The ' // TRIM(msg) // ' field we have does not have the right size' 376 CALL abort_gcm(TRIM(modname), msg, 1) 377 END IF 378 379 END SUBROUTINE get_var_dyn 380 381 !------------------------------------------------------------------------------- 382 383 END SUBROUTINE start_init_dyn 384 385 !------------------------------------------------------------------------------- 386 387 388 !------------------------------------------------------------------------------- 389 390 SUBROUTINE start_inter_3d(var, lon_in, lat_in, lon_in2, lat_in2, pls_in, var3d) 391 392 !------------------------------------------------------------------------------- 393 USE conf_dat_m, ONLY: conf_dat3d 394 USE lmdz_libmath_pch, ONLY: pchsp_95, pchfe_95 395 USE lmdz_libmath, ONLY: minmax 396 397 IMPLICIT NONE 398 !------------------------------------------------------------------------------- 399 ! Arguments: 400 CHARACTER(LEN = *), INTENT(IN) :: var 401 REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) 402 REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) 403 REAL, INTENT(IN) :: pls_in(:, :, :) ! dim (iml,jml,lml) 404 REAL, INTENT(OUT) :: var3d (:, :, :) ! dim (iml,jml,lml) 405 !------------------------------------------------------------------------------- 406 ! Local variables: 407 CHARACTER(LEN = 256) :: modname = 'start_inter_3d' 408 LOGICAL :: skip 409 REAL :: chmin, chmax 410 INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr 411 INTEGER :: n_extrap ! Extrapolated points number 412 REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:) 413 REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:, :, :) 414 REAL, ALLOCATABLE, SAVE :: var_ana3d(:, :, :) 415 !------------------------------------------------------------------------------- 416 iml = assert_eq(SIZE(lon_in), SIZE(lon_in2), SIZE(pls_in, 1), SIZE(var3d, 1), TRIM(modname) // " iml") 417 jml = assert_eq(SIZE(lat_in), SIZE(pls_in, 2), SIZE(var3d, 2), TRIM(modname) // " jml") 418 lml = assert_eq(SIZE(pls_in, 3), SIZE(var3d, 3), TRIM(modname) // " lml"); jml2 = SIZE(lat_in2) 419 420 WRITE(lunout, *)'Going into flinget to extract the 3D field.' 421 IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) 422 CALL flinget(fid_dyn, var, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, var_ana3d) 423 424 !--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS 425 ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) 426 lon_ini(:) = lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini = lon_ini * deg2rad 427 lat_ini(:) = lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini = lat_ini * deg2rad 428 429 !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS 430 ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) 431 CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & 432 lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.) 433 DEALLOCATE(lon_ini, lat_ini) 434 435 !--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro 436 ALLOCATE(var_tmp3d(iml, jml, llm_dyn)) 437 DO il = 1, llm_dyn 438 CALL interp_startvar(var, il==1, lon_rad, lat_rad, var_ana3d(:, :, il), & 439 lon_in, lat_in, lon_in2, lat_in2, var_tmp3d(:, :, il)) 451 440 END DO 452 END DO 453 IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap 454 var3d(iml, :, :) = var3d(1, :, :) 455 456 DO il=1, lml 457 CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax) 458 WRITE(lunout, *)' '//TRIM(var)//' min max l ', il, chmin, chmax 459 END DO 460 461 END SUBROUTINE start_inter_3d 462 463 !------------------------------------------------------------------------------- 464 465 466 !------------------------------------------------------------------------------- 467 468 SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) 469 470 !------------------------------------------------------------------------------- 471 USE inter_barxy_m, ONLY: inter_barxy 472 IMPLICIT NONE 473 !------------------------------------------------------------------------------- 474 ! Arguments: 475 CHARACTER(LEN=*), INTENT(IN) :: nam 476 LOGICAL, INTENT(IN) :: ibeg 477 REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) 478 REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) 479 REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) 480 REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) 481 REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) 482 !------------------------------------------------------------------------------- 483 ! Local variables: 484 CHARACTER(LEN=256) :: modname="interp_startvar" 485 INTEGER :: ii, jj, i1, j1, j2 486 REAL, ALLOCATABLE :: vtmp(:,:) 487 !------------------------------------------------------------------------------- 488 ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") 489 jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") 490 i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") 491 j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1") 492 j2=SIZE(lat2) 493 ALLOCATE(vtmp(i1-1,j1)) 494 IF(ibeg.AND.prt_level>1) THEN 495 WRITE(lunout,*)"---------------------------------------------------------" 496 WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" 497 WRITE(lunout,*)"---------------------------------------------------------" 498 END IF 499 CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) 500 CALL gr_int_dyn(vtmp, varo, i1-1, j1) 501 502 END SUBROUTINE interp_startvar 503 504 !------------------------------------------------------------------------------- 441 DEALLOCATE(lon_rad, lat_rad) 442 443 !--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND 444 ALLOCATE(ax(llm_dyn), ay(llm_dyn), yder(llm_dyn)) 445 ax = lev_dyn(llm_dyn:1:-1) 446 skip = .FALSE. 447 n_extrap = 0 448 DO ij = 1, jml 449 DO ii = 1, iml - 1 450 ay = var_tmp3d(ii, ij, llm_dyn:1:-1) 451 yder = pchsp_95(ax, ay, ibeg = 2, iend = 2, vc_beg = 0., vc_end = 0.) 452 CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & 453 var3d(ii, ij, lml:1:-1), ierr) 454 IF(ierr<0) CALL abort_gcm(TRIM(modname), 'error in pchfe_95', 1) 455 n_extrap = n_extrap + ierr 456 END DO 457 END DO 458 IF(n_extrap/=0) WRITE(lunout, *)TRIM(modname) // " pchfe_95: n_extrap=", n_extrap 459 var3d(iml, :, :) = var3d(1, :, :) 460 461 DO il = 1, lml 462 CALL minmax(iml * jml, var3d(1, 1, il), chmin, chmax) 463 WRITE(lunout, *)' ' // TRIM(var) // ' min max l ', il, chmin, chmax 464 END DO 465 466 END SUBROUTINE start_inter_3d 467 468 !------------------------------------------------------------------------------- 469 470 471 !------------------------------------------------------------------------------- 472 473 SUBROUTINE interp_startvar(nam, ibeg, lon, lat, vari, lon1, lat1, lon2, lat2, varo) 474 475 !------------------------------------------------------------------------------- 476 USE inter_barxy_m, ONLY: inter_barxy 477 IMPLICIT NONE 478 !------------------------------------------------------------------------------- 479 ! Arguments: 480 CHARACTER(LEN = *), INTENT(IN) :: nam 481 LOGICAL, INTENT(IN) :: ibeg 482 REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) 483 REAL, INTENT(IN) :: vari(:, :) ! dim (ii,jj) 484 REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) 485 REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) 486 REAL, INTENT(OUT) :: varo(:, :) ! dim (i1) (j1) 487 !------------------------------------------------------------------------------- 488 ! Local variables: 489 CHARACTER(LEN = 256) :: modname = "interp_startvar" 490 INTEGER :: ii, jj, i1, j1, j2 491 REAL, ALLOCATABLE :: vtmp(:, :) 492 !------------------------------------------------------------------------------- 493 ii = assert_eq(SIZE(lon), SIZE(vari, 1), TRIM(modname) // " ii") 494 jj = assert_eq(SIZE(lat), SIZE(vari, 2), TRIM(modname) // " jj") 495 i1 = assert_eq(SIZE(lon1), SIZE(lon2), SIZE(varo, 1), TRIM(modname) // " i1") 496 j1 = assert_eq(SIZE(lat1), SIZE(varo, 2), TRIM(modname) // " j1") 497 j2 = SIZE(lat2) 498 ALLOCATE(vtmp(i1 - 1, j1)) 499 IF(ibeg.AND.prt_level>1) THEN 500 WRITE(lunout, *)"---------------------------------------------------------" 501 WRITE(lunout, *)"$$$ Interpolation barycentrique pour " // TRIM(nam) // " $$$" 502 WRITE(lunout, *)"---------------------------------------------------------" 503 END IF 504 CALL inter_barxy(lon, lat(:jj - 1), vari, lon2(:i1 - 1), lat2, vtmp) 505 CALL gr_int_dyn(vtmp, varo, i1 - 1, j1) 506 507 END SUBROUTINE interp_startvar 508 509 !------------------------------------------------------------------------------- 505 510 506 511 END MODULE etat0dyn -
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
r5113 r5116 281 281 CALL pbl_surface_init( fder, snsrf, qsurf, tsoil ) 282 282 283 IF (iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1) then283 IF (iflag_pbl>1 .AND. iflag_wake>=1 .AND. iflag_pbl_split >=1) THEN 284 284 delta_tsurf = 0. 285 285 beta_aridity = 0. … … 539 539 do j=2,jmp1-1 540 540 PRINT*,'avant if ',cos(rlatu(j)),coslat0 541 if (cos(rlatu(j))<coslat0) then541 if (cos(rlatu(j))<coslat0) THEN 542 542 ! nb de pts affectes par le filtrage de part et d'autre du pt 543 543 ifiltre=(coslat0/cos(rlatu(j))-1.)/2. … … 548 548 wwf(ifiltre+1)=(coslat0/cos(rlatu(j))-1.)/2.-ifiltre 549 549 do i=1,imp1-1 550 if (masque(i,j)>0.9) then550 if (masque(i,j)>0.9) THEN 551 551 ssz=phis(i,j) 552 552 do ifi=1,ifiltre+1 -
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.f90
r5113 r5116 331 331 USE indice_sol_mod 332 332 USE lmdz_abort_physic, ONLY: abort_physic 333 USE lmdz_libmath, ONLY: minmax 333 334 334 335 IMPLICIT NONE -
LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/test_disvert_m.F90
r5113 r5116 13 13 ! the surface pressure, which sample possible values on Earth. 14 14 15 use exner_hyb_m, only: exner_hyb16 use lmdz_vertical_layers, only: ap,bp,preff17 use comconst_mod, only: kappa, cpp15 use exner_hyb_m, ONLY: exner_hyb 16 use lmdz_vertical_layers, ONLY: ap,bp,preff 17 use comconst_mod, ONLY: kappa, cpp 18 18 USE lmdz_abort_physic, ONLY: abort_physic 19 19 … … 42 42 43 43 ! Are pressure values in the right order? 44 if (any(p(:, :llm) <= p_lay .or. p_lay <= p(:, 2:))) then44 if (any(p(:, :llm) <= p_lay .or. p_lay <= p(:, 2:))) THEN 45 45 ! List details and stop: 46 46 do l = 1, llm 47 47 do i = 1, ngrid 48 if (p(i, l) <= p_lay(i, l)) then48 if (p(i, l) <= p_lay(i, l)) THEN 49 49 print 1000, "ps = ", ps(i) / 100., "hPa, p(level ", l, & 50 50 ") = ", p(i, l) / 100., " hPa <= p(layer ", l, ") = ", & 51 51 p_lay(i, l) / 100., " hPa" 52 52 end if 53 if (p_lay(i, l) <= p(i, l + 1)) then53 if (p_lay(i, l) <= p(i, l + 1)) THEN 54 54 print 1000, "ps = ", ps(i) / 100., "hPa, p(layer ", l, ") = ", & 55 55 p_lay(i, l) / 100., " hPa <= p(level ", l + 1, ") = ", &
Note: See TracChangeset
for help on using the changeset viewer.
