| Line |
Branch |
Exec |
Source |
| 1 |
|
|
MODULE etat0dyn |
| 2 |
|
|
! |
| 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 assert_eq_m, ONLY: assert_eq |
| 38 |
|
|
USE comconst_mod, ONLY: pi, cpp, kappa |
| 39 |
|
|
USE comvert_mod, ONLY: ap, bp, preff, pressure_exner |
| 40 |
|
|
USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itau_phy, start_time |
| 41 |
|
|
|
| 42 |
|
|
IMPLICIT NONE |
| 43 |
|
|
|
| 44 |
|
|
PRIVATE |
| 45 |
|
|
PUBLIC :: etat0dyn_netcdf |
| 46 |
|
|
|
| 47 |
|
|
include "iniprint.h" |
| 48 |
|
|
include "dimensions.h" |
| 49 |
|
|
include "paramet.h" |
| 50 |
|
|
include "comgeom2.h" |
| 51 |
|
|
include "comdissnew.h" |
| 52 |
|
|
REAL, SAVE :: deg2rad |
| 53 |
|
|
INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn |
| 54 |
|
|
REAL, ALLOCATABLE, SAVE :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:) |
| 55 |
|
|
CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc' |
| 56 |
|
|
|
| 57 |
|
|
CONTAINS |
| 58 |
|
|
|
| 59 |
|
|
!------------------------------------------------------------------------------- |
| 60 |
|
|
! |
| 61 |
|
✗ |
SUBROUTINE etat0dyn_netcdf(masque, phis) |
| 62 |
|
|
! |
| 63 |
|
|
!------------------------------------------------------------------------------- |
| 64 |
|
|
! Purpose: Create dynamical initial states. |
| 65 |
|
|
!------------------------------------------------------------------------------- |
| 66 |
|
|
! Notes: 1) This routine is designed to work for Earth |
| 67 |
|
|
! 2) If masque(:,:)/=-99999., masque and phis are already known. |
| 68 |
|
|
! Otherwise: compute it. |
| 69 |
|
|
!------------------------------------------------------------------------------- |
| 70 |
|
|
USE control_mod |
| 71 |
|
|
USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz |
| 72 |
|
|
USE regr_pr_o3_m, ONLY: regr_pr_o3 |
| 73 |
|
|
USE press_coefoz_m, ONLY: press_coefoz |
| 74 |
|
|
USE exner_hyb_m, ONLY: exner_hyb |
| 75 |
|
|
USE exner_milieu_m, ONLY: exner_milieu |
| 76 |
|
|
USE infotrac, ONLY: nqtot, tname |
| 77 |
|
|
USE filtreg_mod |
| 78 |
|
|
IMPLICIT NONE |
| 79 |
|
|
!------------------------------------------------------------------------------- |
| 80 |
|
|
! Arguments: |
| 81 |
|
|
REAL, INTENT(INOUT) :: masque(iip1,jjp1) !--- Land-ocean mask |
| 82 |
|
|
REAL, INTENT(INOUT) :: phis (iip1,jjp1) !--- Ground geopotential |
| 83 |
|
|
!------------------------------------------------------------------------------- |
| 84 |
|
|
! Local variables: |
| 85 |
|
|
CHARACTER(LEN=256) :: modname, fmt |
| 86 |
|
|
INTEGER :: i, j, l, ji, itau, iday |
| 87 |
|
|
REAL :: xpn, xps, time, phystep |
| 88 |
|
|
REAL, DIMENSION(iip1,jjp1) :: psol |
| 89 |
|
|
REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d |
| 90 |
|
|
REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd |
| 91 |
|
|
REAL, DIMENSION(iip1,jjp1,llm) :: pk, pls, y, masse |
| 92 |
|
|
REAL, DIMENSION(iip1,jjm ,llm) :: vvent |
| 93 |
|
|
REAL, DIMENSION(ip1jm ,llm) :: pbarv |
| 94 |
|
|
REAL, DIMENSION(ip1jmp1 ,llm) :: pbaru, phi, w |
| 95 |
|
|
REAL, DIMENSION(ip1jmp1) :: pks |
| 96 |
|
|
REAL, DIMENSION(iim) :: xppn, xpps |
| 97 |
|
|
REAL, ALLOCATABLE :: q3d(:,:,:,:) |
| 98 |
|
|
!------------------------------------------------------------------------------- |
| 99 |
|
✗ |
modname='etat0dyn_netcdf' |
| 100 |
|
|
|
| 101 |
|
✗ |
deg2rad = pi/180.0 |
| 102 |
|
✗ |
y(:,:,:)=0 !ym warning unitialized variable |
| 103 |
|
|
|
| 104 |
|
|
! Compute psol AND tsol, knowing phis. |
| 105 |
|
|
!******************************************************************************* |
| 106 |
|
✗ |
CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol) |
| 107 |
|
|
|
| 108 |
|
|
! Mid-levels pressure computation |
| 109 |
|
|
!******************************************************************************* |
| 110 |
|
✗ |
CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d |
| 111 |
|
✗ |
IF(pressure_exner) THEN !--- Update pk, pks |
| 112 |
|
✗ |
CALL exner_hyb (ip1jmp1,psol,p3d,pks,pk) |
| 113 |
|
|
ELSE |
| 114 |
|
✗ |
CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk) |
| 115 |
|
|
END IF |
| 116 |
|
✗ |
pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa) !--- Update pls |
| 117 |
|
|
|
| 118 |
|
|
! Update uvent, vvent, t3d and tpot |
| 119 |
|
|
!******************************************************************************* |
| 120 |
|
✗ |
uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0 |
| 121 |
|
✗ |
CALL startget_dyn3d('u' ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv) |
| 122 |
|
|
CALL startget_dyn3d('v' ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent, & |
| 123 |
|
✗ |
& rlonu,rlatu(:jjm)) |
| 124 |
|
✗ |
CALL startget_dyn3d('t' ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv) |
| 125 |
|
✗ |
tpot(:,:,:)=t3d(:,:,:) |
| 126 |
|
✗ |
CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv) |
| 127 |
|
|
|
| 128 |
|
✗ |
WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:)) |
| 129 |
|
✗ |
WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:)) |
| 130 |
|
|
|
| 131 |
|
|
! Humidity at saturation computation |
| 132 |
|
|
!******************************************************************************* |
| 133 |
|
✗ |
WRITE(lunout,*) 'avant q_sat' |
| 134 |
|
✗ |
CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat) |
| 135 |
|
✗ |
WRITE(lunout,*) 'apres q_sat' |
| 136 |
|
✗ |
WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:)) |
| 137 |
|
|
! WRITE(lunout,*) 'QSAT :',qsat(10,20,:) |
| 138 |
|
✗ |
qd (:,:,:) = 0.0 |
| 139 |
|
✗ |
CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv) |
| 140 |
|
✗ |
ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:) |
| 141 |
|
✗ |
CALL flinclo(fid_dyn) |
| 142 |
|
|
|
| 143 |
|
|
! Parameterization of ozone chemistry: |
| 144 |
|
|
!******************************************************************************* |
| 145 |
|
|
! Look for ozone tracer: |
| 146 |
|
✗ |
DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO |
| 147 |
|
✗ |
IF(i/=nqtot+1) THEN |
| 148 |
|
✗ |
CALL regr_lat_time_coefoz |
| 149 |
|
✗ |
CALL press_coefoz |
| 150 |
|
✗ |
CALL regr_pr_o3(p3d, q3d(:,:,:,i)) |
| 151 |
|
✗ |
q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29. !--- Mole->mass fraction |
| 152 |
|
|
END IF |
| 153 |
|
✗ |
q3d(iip1,:,:,:)=q3d(1,:,:,:) |
| 154 |
|
|
|
| 155 |
|
|
! Writing |
| 156 |
|
|
!******************************************************************************* |
| 157 |
|
|
CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, & |
| 158 |
|
✗ |
tetatemp, vert_prof_dissip) |
| 159 |
|
✗ |
WRITE(lunout,*)'sortie inidissip' |
| 160 |
|
✗ |
itau=0 |
| 161 |
|
✗ |
itau_dyn=0 |
| 162 |
|
✗ |
itau_phy=0 |
| 163 |
|
✗ |
iday=dayref+itau/day_step |
| 164 |
|
✗ |
time=FLOAT(itau-(iday-dayref)*day_step)/day_step |
| 165 |
|
✗ |
IF(time>1.) THEN |
| 166 |
|
✗ |
time=time-1 |
| 167 |
|
✗ |
iday=iday+1 |
| 168 |
|
|
END IF |
| 169 |
|
✗ |
day_ref=dayref |
| 170 |
|
✗ |
annee_ref=anneeref |
| 171 |
|
✗ |
CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi ) |
| 172 |
|
✗ |
WRITE(lunout,*)'sortie geopot' |
| 173 |
|
|
CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, & |
| 174 |
|
✗ |
phi, w, pbaru, pbarv, time+iday-dayref) |
| 175 |
|
✗ |
WRITE(lunout,*)'sortie caldyn0' |
| 176 |
|
✗ |
start_time = 0. |
| 177 |
|
✗ |
CALL dynredem0( "start.nc", dayref, phis) |
| 178 |
|
✗ |
WRITE(lunout,*)'sortie dynredem0' |
| 179 |
|
✗ |
CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol) |
| 180 |
|
✗ |
WRITE(lunout,*)'sortie dynredem1' |
| 181 |
|
✗ |
CALL histclo() |
| 182 |
|
|
|
| 183 |
|
✗ |
END SUBROUTINE etat0dyn_netcdf |
| 184 |
|
|
! |
| 185 |
|
|
!------------------------------------------------------------------------------- |
| 186 |
|
|
|
| 187 |
|
|
|
| 188 |
|
|
!------------------------------------------------------------------------------- |
| 189 |
|
|
! |
| 190 |
|
✗ |
SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar,& |
| 191 |
|
✗ |
champ, lon_in2, lat_in2) |
| 192 |
|
|
!------------------------------------------------------------------------------- |
| 193 |
|
|
IMPLICIT NONE |
| 194 |
|
|
!=============================================================================== |
| 195 |
|
|
! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R |
| 196 |
|
|
! (3D fields) of file dynfname. |
| 197 |
|
|
!------------------------------------------------------------------------------- |
| 198 |
|
|
! Note: An input auxilliary field "workvar" has to be specified in two cases: |
| 199 |
|
|
! * for "q": the saturated humidity. |
| 200 |
|
|
! * for "tpot": the Exner function. |
| 201 |
|
|
!=============================================================================== |
| 202 |
|
|
! Arguments: |
| 203 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: var |
| 204 |
|
|
REAL, INTENT(IN) :: lon_in(:) ! dim (iml) |
| 205 |
|
|
REAL, INTENT(IN) :: lat_in(:) ! dim (jml) |
| 206 |
|
|
REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml) |
| 207 |
|
|
REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml) |
| 208 |
|
|
REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml) |
| 209 |
|
|
REAL, INTENT(IN) :: lon_in2(:) ! dim (iml) |
| 210 |
|
|
REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2) |
| 211 |
|
|
!------------------------------------------------------------------------------- |
| 212 |
|
|
! Local variables: |
| 213 |
|
|
CHARACTER(LEN=10) :: vname |
| 214 |
|
|
CHARACTER(LEN=256) :: msg, modname="startget_dyn3d" |
| 215 |
|
|
INTEGER :: iml, jml, jml2, lml, il |
| 216 |
|
|
REAL :: xppn, xpps |
| 217 |
|
|
!------------------------------------------------------------------------------- |
| 218 |
|
|
iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1), & |
| 219 |
|
✗ |
& SIZE(lon_in2)], TRIM(modname)//" iml") |
| 220 |
|
|
jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2), & |
| 221 |
|
✗ |
& TRIM(modname)//" jml") |
| 222 |
|
|
lml=assert_eq( SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3), & |
| 223 |
|
✗ |
& TRIM(modname)//" lml") |
| 224 |
|
|
jml2=SIZE(lat_in2) |
| 225 |
|
|
|
| 226 |
|
|
!--- CHECK IF THE FIELD IS KNOWN |
| 227 |
|
✗ |
SELECT CASE(var) |
| 228 |
|
✗ |
CASE('u'); vname='U' |
| 229 |
|
✗ |
CASE('v'); vname='V' |
| 230 |
|
✗ |
CASE('t'); vname='TEMP' |
| 231 |
|
✗ |
CASE('q'); vname='R'; msg='humidity as the saturated humidity' |
| 232 |
|
✗ |
CASE('tpot'); msg='potential temperature as the Exner function' |
| 233 |
|
✗ |
CASE DEFAULT; msg='No rule to extract variable '//TRIM(var) |
| 234 |
|
✗ |
CALL abort_gcm(modname,TRIM(msg)//' from any data set',1) |
| 235 |
|
|
END SELECT |
| 236 |
|
|
|
| 237 |
|
|
!--- CHECK IF SOMETHING IS MISSING |
| 238 |
|
✗ |
IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN |
| 239 |
|
✗ |
msg='Could not compute '//TRIM(msg)//' is missing or constant.' |
| 240 |
|
✗ |
CALL abort_gcm(modname,TRIM(msg),1) |
| 241 |
|
|
END IF |
| 242 |
|
|
|
| 243 |
|
|
!--- INTERPOLATE 3D FIELD IF NEEDED |
| 244 |
|
✗ |
IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2, & |
| 245 |
|
✗ |
lat_in2,pls,champ) |
| 246 |
|
|
|
| 247 |
|
|
!--- COMPUTE THE REQUIRED FILED |
| 248 |
|
|
SELECT CASE(var) |
| 249 |
|
✗ |
CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO |
| 250 |
|
✗ |
champ(iml,:,:)=champ(1,:,:) !--- Eastward wind |
| 251 |
|
|
|
| 252 |
|
✗ |
CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO |
| 253 |
|
✗ |
champ(iml,:,:)=champ(1,:,:) !--- Northward wind |
| 254 |
|
|
|
| 255 |
|
|
CASE('tpot','q') |
| 256 |
|
✗ |
IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature |
| 257 |
|
✗ |
ELSE; champ=champ*.01*workvar !--- Relative humidity |
| 258 |
|
✗ |
WHERE(champ<0.) champ=1.0E-10 |
| 259 |
|
|
END IF |
| 260 |
|
✗ |
DO il=1,lml |
| 261 |
|
✗ |
xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln |
| 262 |
|
✗ |
xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols |
| 263 |
|
✗ |
champ(:,1 ,il) = xppn |
| 264 |
|
✗ |
champ(:,jml,il) = xpps |
| 265 |
|
|
END DO |
| 266 |
|
|
END SELECT |
| 267 |
|
|
|
| 268 |
|
✗ |
END SUBROUTINE startget_dyn3d |
| 269 |
|
|
! |
| 270 |
|
|
!------------------------------------------------------------------------------- |
| 271 |
|
|
|
| 272 |
|
|
|
| 273 |
|
|
!------------------------------------------------------------------------------- |
| 274 |
|
|
! |
| 275 |
|
✗ |
SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol) |
| 276 |
|
|
! |
| 277 |
|
|
!------------------------------------------------------------------------------- |
| 278 |
|
|
IMPLICIT NONE |
| 279 |
|
|
!=============================================================================== |
| 280 |
|
|
! Purpose: Compute psol, knowing phis. |
| 281 |
|
|
!=============================================================================== |
| 282 |
|
|
! Arguments: |
| 283 |
|
|
REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml) |
| 284 |
|
|
REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) |
| 285 |
|
|
REAL, INTENT(IN) :: zs (:,:) ! dim (iml,jml) |
| 286 |
|
|
REAL, INTENT(OUT) :: psol(:,:) ! dim (iml,jml) |
| 287 |
|
|
!------------------------------------------------------------------------------- |
| 288 |
|
|
! Local variables: |
| 289 |
|
|
CHARACTER(LEN=256) :: modname='start_init_dyn' |
| 290 |
|
|
REAL :: date, dt |
| 291 |
|
|
INTEGER :: iml, jml, jml2, itau(1) |
| 292 |
|
✗ |
REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:,:) |
| 293 |
|
✗ |
REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:) |
| 294 |
|
✗ |
REAL, ALLOCATABLE :: z(:,:), ps(:,:), ts(:,:) |
| 295 |
|
|
!------------------------------------------------------------------------------- |
| 296 |
|
|
iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2), & |
| 297 |
|
✗ |
& TRIM(modname)//" iml") |
| 298 |
|
✗ |
jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml") |
| 299 |
|
|
jml2=SIZE(lat_in2) |
| 300 |
|
|
|
| 301 |
|
✗ |
WRITE(lunout,*) 'Opening the surface analysis' |
| 302 |
|
✗ |
CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn) |
| 303 |
|
✗ |
WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn |
| 304 |
|
|
|
| 305 |
|
✗ |
ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn)) |
| 306 |
|
✗ |
ALLOCATE(levdyn_ini(llm_dyn)) |
| 307 |
|
|
CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, & |
| 308 |
|
✗ |
lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn) |
| 309 |
|
|
|
| 310 |
|
|
!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS |
| 311 |
|
✗ |
ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn)) |
| 312 |
|
✗ |
lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad |
| 313 |
|
✗ |
lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad |
| 314 |
|
|
|
| 315 |
|
✗ |
ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn)) |
| 316 |
|
✗ |
CALL get_var_dyn('Z',z) !--- SURFACE GEOPOTENTIAL |
| 317 |
|
✗ |
CALL get_var_dyn('SP',ps) !--- SURFACE PRESSURE |
| 318 |
|
✗ |
CALL get_var_dyn('ST',ts) !--- SURFACE TEMPERATURE |
| 319 |
|
|
! CALL flinclo(fid_dyn) |
| 320 |
|
✗ |
DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini) |
| 321 |
|
|
|
| 322 |
|
|
!--- PSOL IS COMPUTED IN PASCALS |
| 323 |
|
|
psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0 & |
| 324 |
|
✗ |
& /ts(:iml-1,:)) |
| 325 |
|
✗ |
psol(iml,:)=psol(1,:) |
| 326 |
|
✗ |
DEALLOCATE(z,ps,ts) |
| 327 |
|
✗ |
psol(:,1 )=SUM(aire(1:iml-1,1 )*psol(1:iml-1,1 ))/apoln !--- NORTH POLE |
| 328 |
|
✗ |
psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols !--- SOUTH POLE |
| 329 |
|
|
|
| 330 |
|
|
CONTAINS |
| 331 |
|
|
|
| 332 |
|
|
!------------------------------------------------------------------------------- |
| 333 |
|
|
! |
| 334 |
|
✗ |
SUBROUTINE get_var_dyn(title,field) |
| 335 |
|
|
! |
| 336 |
|
|
!------------------------------------------------------------------------------- |
| 337 |
|
|
USE conf_dat_m, ONLY: conf_dat2d |
| 338 |
|
|
IMPLICIT NONE |
| 339 |
|
|
!------------------------------------------------------------------------------- |
| 340 |
|
|
! Arguments: |
| 341 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: title |
| 342 |
|
|
REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:) |
| 343 |
|
|
!------------------------------------------------------------------------------- |
| 344 |
|
|
! Local variables: |
| 345 |
|
|
CHARACTER(LEN=256) :: msg |
| 346 |
|
|
INTEGER :: tllm |
| 347 |
|
|
!------------------------------------------------------------------------------- |
| 348 |
|
✗ |
SELECT CASE(title) |
| 349 |
|
✗ |
CASE('Z'); tllm=0; msg='geopotential' |
| 350 |
|
✗ |
CASE('SP'); tllm=0; msg='surface pressure' |
| 351 |
|
✗ |
CASE('ST'); tllm=llm_dyn; msg='temperature' |
| 352 |
|
|
END SELECT |
| 353 |
|
✗ |
IF(.NOT.ALLOCATED(field)) THEN |
| 354 |
|
✗ |
ALLOCATE(field(iml,jml)) |
| 355 |
|
✗ |
CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana) |
| 356 |
|
✗ |
CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.) |
| 357 |
|
|
CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana, & |
| 358 |
|
✗ |
lon_in, lat_in, lon_in2, lat_in2, field) |
| 359 |
|
✗ |
ELSE IF(SIZE(field)/=SIZE(z)) THEN |
| 360 |
|
✗ |
msg='The '//TRIM(msg)//' field we have does not have the right size' |
| 361 |
|
✗ |
CALL abort_gcm(TRIM(modname),msg,1) |
| 362 |
|
|
END IF |
| 363 |
|
|
|
| 364 |
|
✗ |
END SUBROUTINE get_var_dyn |
| 365 |
|
|
! |
| 366 |
|
|
!------------------------------------------------------------------------------- |
| 367 |
|
|
|
| 368 |
|
|
END SUBROUTINE start_init_dyn |
| 369 |
|
|
! |
| 370 |
|
|
!------------------------------------------------------------------------------- |
| 371 |
|
|
|
| 372 |
|
|
|
| 373 |
|
|
!------------------------------------------------------------------------------- |
| 374 |
|
|
! |
| 375 |
|
✗ |
SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d) |
| 376 |
|
|
! |
| 377 |
|
|
!------------------------------------------------------------------------------- |
| 378 |
|
|
USE conf_dat_m, ONLY: conf_dat3d |
| 379 |
|
|
USE pchsp_95_m, ONLY: pchsp_95 |
| 380 |
|
|
USE pchfe_95_m, ONLY: pchfe_95 |
| 381 |
|
|
IMPLICIT NONE |
| 382 |
|
|
!------------------------------------------------------------------------------- |
| 383 |
|
|
! Arguments: |
| 384 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: var |
| 385 |
|
|
REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml) |
| 386 |
|
|
REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2) |
| 387 |
|
|
REAL, INTENT(IN) :: pls_in(:,:,:) ! dim (iml,jml,lml) |
| 388 |
|
|
REAL, INTENT(OUT) :: var3d (:,:,:) ! dim (iml,jml,lml) |
| 389 |
|
|
!------------------------------------------------------------------------------- |
| 390 |
|
|
! Local variables: |
| 391 |
|
|
CHARACTER(LEN=256) :: modname='start_inter_3d' |
| 392 |
|
|
LOGICAL :: skip |
| 393 |
|
|
REAL :: chmin, chmax |
| 394 |
|
|
INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr |
| 395 |
|
|
INTEGER :: n_extrap ! Extrapolated points number |
| 396 |
|
✗ |
REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:) |
| 397 |
|
✗ |
REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:) |
| 398 |
|
|
REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:) |
| 399 |
|
|
!------------------------------------------------------------------------------- |
| 400 |
|
✗ |
iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml") |
| 401 |
|
✗ |
jml=assert_eq(SIZE(lat_in), SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml") |
| 402 |
|
✗ |
lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2) |
| 403 |
|
|
|
| 404 |
|
✗ |
WRITE(lunout, *)'Going into flinget to extract the 3D field.' |
| 405 |
|
✗ |
IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn)) |
| 406 |
|
✗ |
CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d) |
| 407 |
|
|
|
| 408 |
|
|
!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS |
| 409 |
|
✗ |
ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn)) |
| 410 |
|
✗ |
lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad |
| 411 |
|
✗ |
lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad |
| 412 |
|
|
|
| 413 |
|
|
!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS |
| 414 |
|
✗ |
ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn)) |
| 415 |
|
|
CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, & |
| 416 |
|
✗ |
lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.) |
| 417 |
|
✗ |
DEALLOCATE(lon_ini, lat_ini) |
| 418 |
|
|
|
| 419 |
|
|
!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro |
| 420 |
|
✗ |
ALLOCATE(var_tmp3d(iml,jml,llm_dyn)) |
| 421 |
|
✗ |
DO il = 1,llm_dyn |
| 422 |
|
|
CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il), & |
| 423 |
|
✗ |
lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il)) |
| 424 |
|
|
END DO |
| 425 |
|
✗ |
DEALLOCATE(lon_rad, lat_rad) |
| 426 |
|
|
|
| 427 |
|
|
!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND |
| 428 |
|
✗ |
ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn)) |
| 429 |
|
✗ |
ax = lev_dyn(llm_dyn:1:-1) |
| 430 |
|
✗ |
skip = .FALSE. |
| 431 |
|
✗ |
n_extrap = 0 |
| 432 |
|
✗ |
DO ij=1, jml |
| 433 |
|
✗ |
DO ii=1, iml-1 |
| 434 |
|
✗ |
ay = var_tmp3d(ii, ij, llm_dyn:1:-1) |
| 435 |
|
✗ |
yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.) |
| 436 |
|
|
CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), & |
| 437 |
|
✗ |
var3d(ii, ij, lml:1:-1), ierr) |
| 438 |
|
✗ |
IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1) |
| 439 |
|
✗ |
n_extrap = n_extrap + ierr |
| 440 |
|
|
END DO |
| 441 |
|
|
END DO |
| 442 |
|
✗ |
IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap |
| 443 |
|
✗ |
var3d(iml, :, :) = var3d(1, :, :) |
| 444 |
|
|
|
| 445 |
|
✗ |
DO il=1, lml |
| 446 |
|
✗ |
CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax) |
| 447 |
|
✗ |
WRITE(lunout, *)' '//TRIM(var)//' min max l ', il, chmin, chmax |
| 448 |
|
|
END DO |
| 449 |
|
|
|
| 450 |
|
✗ |
END SUBROUTINE start_inter_3d |
| 451 |
|
|
! |
| 452 |
|
|
!------------------------------------------------------------------------------- |
| 453 |
|
|
|
| 454 |
|
|
|
| 455 |
|
|
!------------------------------------------------------------------------------- |
| 456 |
|
|
! |
| 457 |
|
✗ |
SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo) |
| 458 |
|
|
! |
| 459 |
|
|
!------------------------------------------------------------------------------- |
| 460 |
|
|
USE inter_barxy_m, ONLY: inter_barxy |
| 461 |
|
|
IMPLICIT NONE |
| 462 |
|
|
!------------------------------------------------------------------------------- |
| 463 |
|
|
! Arguments: |
| 464 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: nam |
| 465 |
|
|
LOGICAL, INTENT(IN) :: ibeg |
| 466 |
|
|
REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj) |
| 467 |
|
|
REAL, INTENT(IN) :: vari(:,:) ! dim (ii,jj) |
| 468 |
|
|
REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1) |
| 469 |
|
|
REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2) |
| 470 |
|
|
REAL, INTENT(OUT) :: varo(:,:) ! dim (i1) (j1) |
| 471 |
|
|
!------------------------------------------------------------------------------- |
| 472 |
|
|
! Local variables: |
| 473 |
|
|
CHARACTER(LEN=256) :: modname="interp_startvar" |
| 474 |
|
|
INTEGER :: ii, jj, i1, j1, j2 |
| 475 |
|
✗ |
REAL, ALLOCATABLE :: vtmp(:,:) |
| 476 |
|
|
!------------------------------------------------------------------------------- |
| 477 |
|
✗ |
ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii") |
| 478 |
|
✗ |
jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj") |
| 479 |
|
✗ |
i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1") |
| 480 |
|
✗ |
j1=assert_eq(SIZE(lat1), SIZE(varo,2),TRIM(modname)//" j1") |
| 481 |
|
|
j2=SIZE(lat2) |
| 482 |
|
✗ |
ALLOCATE(vtmp(i1-1,j1)) |
| 483 |
|
✗ |
IF(ibeg.AND.prt_level>1) THEN |
| 484 |
|
✗ |
WRITE(lunout,*)"---------------------------------------------------------" |
| 485 |
|
✗ |
WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$" |
| 486 |
|
✗ |
WRITE(lunout,*)"---------------------------------------------------------" |
| 487 |
|
|
END IF |
| 488 |
|
✗ |
CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp) |
| 489 |
|
✗ |
CALL gr_int_dyn(vtmp, varo, i1-1, j1) |
| 490 |
|
|
|
| 491 |
|
✗ |
END SUBROUTINE interp_startvar |
| 492 |
|
|
! |
| 493 |
|
|
!------------------------------------------------------------------------------- |
| 494 |
|
|
|
| 495 |
|
|
END MODULE etat0dyn |
| 496 |
|
|
! |
| 497 |
|
|
!******************************************************************************* |
| 498 |
|
|
|