source: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90 @ 2333

Last change on this file since 2333 was 2331, checked in by lguez, 9 years ago

Fixed regression from revision 2315: comvert.h was replaced by
vertical_layers_mod in test_disvert, but variables ap, bp, preff of
vertical_layers_mod were not defined. So, in main program ce0l, moved
call to test_disvert after call to Init_Phys_lmdz, and inserted in
between them calls to infotrac_init and iniphysiq (required). Had then
to remove the call to infotrac_init in etat0dyn_netcdf. In main
program ce0l, had to remove the call to InitComgeomphy? since this is
done in iniphysiq.

In main program ce0l: no need to use indice_sol_mod; removed
preprocessor tests on CPP_PHYS in ce0l.
File size: 32.1 KB
RevLine 
[2293]1MODULE 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,&
[2299]14!                          champ, lon_in2, lat_in2, ibar)
[2293]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#ifdef CPP_PHYS
39  USE indice_sol_mod, ONLY: epsfra
40#endif
41  IMPLICIT NONE
42
43  PRIVATE
44  PUBLIC :: etat0dyn_netcdf
45
46  include "iniprint.h"
47  include "dimensions.h"
48  include "paramet.h"
49  include "comgeom2.h"
50  include "comvert.h"
51  include "comconst.h"
52  include "temps.h"
53  include "comdissnew.h"
54  include "serre.h"
55  REAL, SAVE :: deg2rad
56#ifndef CPP_PHYS
57  REAL, SAVE :: epsfra= 1.E-5
58#endif
59  INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
60  REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
61  CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
62  CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'
63
64CONTAINS
65
66!-------------------------------------------------------------------------------
67!
68SUBROUTINE etat0dyn_netcdf(ib, masque, phis)
69!
70!-------------------------------------------------------------------------------
71! Purpose: Create dynamical initial states.
72!-------------------------------------------------------------------------------
73! Notes:  1) This routine is designed to work for Earth
74!         2) If masque(:,:)/=-99999., masque and phis are already known.
75!         Otherwise: read it from ocean model file (coupled run) or compute it.
76!-------------------------------------------------------------------------------
77  USE control_mod
78#ifdef CPP_PHYS
79  USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
80  USE regr_pr_o3_m,   ONLY: regr_pr_o3
81  USE press_coefoz_m, ONLY: press_coefoz
82#endif
83  USE exner_hyb_m,    ONLY: exner_hyb
84  USE exner_milieu_m, ONLY: exner_milieu
[2331]85  USE infotrac, only: NQTOT, TNAME
[2293]86  USE filtreg_mod
87  IMPLICIT NONE
88!-------------------------------------------------------------------------------
89! Arguments:
90  LOGICAL, INTENT(IN)    :: ib                  !--- Barycentric interpolation
91  REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
92  REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
93!-------------------------------------------------------------------------------
94! Local variables:
95  CHARACTER(LEN=256) :: modname, fmt
96  INTEGER            :: i, j, l, ji, itau, iday
97  REAL               :: xpn, xps, time, phystep
98  REAL, DIMENSION(iip1,jjp1)       :: psol, masque_tmp
99  REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
100  REAL, DIMENSION(iip1,jjp1,llm)   :: uvent, t3d, tpot, qsat, qd
101  REAL, DIMENSION(iip1,jjp1,llm)   :: pk, pls, y, masse
102  REAL, DIMENSION(iip1,jjm ,llm)   :: vvent
103  REAL, DIMENSION(ip1jm    ,llm)   :: pbarv
104  REAL, DIMENSION(ip1jmp1  ,llm)   :: pbaru, phi, w
105  REAL, DIMENSION(ip1jmp1)         :: pks
106  REAL, DIMENSION(iim)             :: xppn, xpps
107  REAL, ALLOCATABLE                :: q3d(:,:,:,:)
108!-------------------------------------------------------------------------------
109  modname='etat0dyn_netcdf'
110
111  deg2rad = pi/180.0
112
113! Initializations for tracers and filter
114!*******************************************************************************
115  CALL inifilr()
116
117! Compute ground geopotential and possibly the mask.
118!*******************************************************************************
119  masque_tmp(:,:)=masque(:,:)
120  CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)
121  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
122  IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
123    masque=masque_tmp
124    IF(prt_level>=1) THEN
125      WRITE(lunout,*)'BUILT MASK :'
126      WRITE(lunout,fmt) NINT(masque)
127    END IF
128    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
129    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
130  END IF
131
132! Compute psol AND tsol, knowing phis.
133!*******************************************************************************
134  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib, phis, psol)
135
136! Mid-levels pressure computation
137!*******************************************************************************
138  CALL pression(ip1jmp1, ap, bp, psol, p3d)             !--- Update p3d
139  IF(pressure_exner) THEN                               !--- Update pk, pks
140    CALL exner_hyb   (ip1jmp1,psol,p3d,pks,pk)
141  ELSE
142    CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
143  END IF
144  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)          !--- Update pls
145
146! Update uvent, vvent, t3d and tpot
147!*******************************************************************************
[2299]148  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
149  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
150  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
[2293]151 &                           rlonu,rlatu(:jjm),ib)
[2299]152  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
153  tpot(:,:,:)=t3d(:,:,:)
154  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
[2293]155
156  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
157  WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
158
159! Humidity at saturation computation
160!*******************************************************************************
161  WRITE(lunout,*) 'avant q_sat'
162  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
163  WRITE(lunout,*) 'apres q_sat'
164  WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
165!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
166  qd (:,:,:) = 0.0
[2299]167  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
168  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
[2293]169  CALL flinclo(fid_dyn)
170
171#ifdef CPP_PHYS
172! Parameterization of ozone chemistry:
173!*******************************************************************************
174! Look for ozone tracer:
175  DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
176  IF(i/=nqtot+1) THEN
177    CALL regr_lat_time_coefoz
178    CALL press_coefoz
179    CALL regr_pr_o3(p3d, q3d(:,:,:,i))
180    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
181  END IF
[2299]182
[2293]183#endif
[2299]184  q3d(iip1,:,:,:)=q3d(1,:,:,:)
[2293]185
186! Intermediate computation
187!*******************************************************************************
188  CALL massdair(p3d,masse)
189  WRITE(lunout,*)' ALPHAX ',alphax
190  DO l=1,llm
191    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
192    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
193    xpn=SUM(xppn)/apoln
194    xps=SUM(xpps)/apols
195    masse(:,1   ,l)=xpn
196    masse(:,jjp1,l)=xps
197  END DO
198
199! Writing
200!*******************************************************************************
201  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot,   &
202       tetatemp, vert_prof_dissip)
203  WRITE(lunout,*)'sortie inidissip'
204  itau=0
205  itau_dyn=0
206  itau_phy=0
207  iday=dayref+itau/day_step
208  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
209  IF(time>1.) THEN
210   time=time-1
211   iday=iday+1
212  END IF
213  day_ref=dayref
214  annee_ref=anneeref
215  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
216  WRITE(lunout,*)'sortie geopot'
217  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
218                phi,  w, pbaru, pbarv, time+iday-dayref)
219  WRITE(lunout,*)'sortie caldyn0'     
220  CALL dynredem0( "start.nc", dayref, phis)
221  WRITE(lunout,*)'sortie dynredem0'
222  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
223  WRITE(lunout,*)'sortie dynredem1'
224  CALL histclo()
225
226END SUBROUTINE etat0dyn_netcdf
227!
228!-------------------------------------------------------------------------------
229
230
231!-------------------------------------------------------------------------------
232!
[2299]233SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
234                        champ, lon_in2, lat_in2, ibar)
[2293]235!-------------------------------------------------------------------------------
236  IMPLICIT NONE
237!===============================================================================
238! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
239!     (3D fields) of file dynfname.
240!-------------------------------------------------------------------------------
241! Note: An input auxilliary field "workvar" has to be specified in two cases:
[2299]242!     * for "q":    the saturated humidity.
243!     * for "tpot": the Exner function.
[2293]244!===============================================================================
245! Arguments:
246  CHARACTER(LEN=*), INTENT(IN)    :: var
247  REAL,             INTENT(IN)    :: lon_in(:)        ! dim (iml)
248  REAL,             INTENT(IN)    :: lat_in(:)        ! dim (jml)
249  REAL,             INTENT(IN)    :: pls    (:, :, :) ! dim (iml, jml, lml)
250  REAL,             INTENT(IN)    :: workvar(:, :, :) ! dim (iml, jml, lml)
251  REAL,             INTENT(INOUT) :: champ  (:, :, :) ! dim (iml, jml, lml)
252  REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
253  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
254  LOGICAL,          INTENT(IN)    :: ibar
255!-------------------------------------------------------------------------------
256! Local variables:
257  CHARACTER(LEN=10)  :: vname
258  CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
259  INTEGER            :: iml, jml, jml2, lml, il
260  REAL               :: xppn, xpps
261!-------------------------------------------------------------------------------
[2299]262  iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),       &
263     &                                    SIZE(lon_in2)], TRIM(modname)//" iml")
264  jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),       &
265     &                                                    TRIM(modname)//" jml")
266  lml=assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),       &
267     &                                                    TRIM(modname)//" lml")
268  jml2=SIZE(lat_in2)
[2293]269
[2299]270!--- CHECK IF THE FIELD IS KNOWN
271   SELECT CASE(var)
272    CASE('u');    vname='U'
273    CASE('v');    vname='V'
274    CASE('t');    vname='TEMP'
275    CASE('q');    vname='R';    msg='humidity as the saturated humidity'
276    CASE('tpot'); msg='potential temperature as the Exner function'
277    CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
278      CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
279  END SELECT
[2293]280
[2299]281!--- CHECK IF SOMETHING IS MISSING
282  IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
283    msg='Could not compute '//TRIM(msg)//' is missing or constant.'
284    CALL abort_gcm(modname,TRIM(msg),1)
285  END IF
[2293]286
[2299]287!--- INTERPOLATE 3D FIELD IF NEEDED
288  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
289                                                  lat_in2,pls,champ,ibar)
[2293]290
[2299]291!--- COMPUTE THE REQUIRED FILED
292  SELECT CASE(var)
293    CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
294      champ(iml,:,:)=champ(1,:,:)                   !--- Eastward wind
[2293]295
[2299]296    CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
297      champ(iml,:,:)=champ(1,:,:)                   !--- Northward wind
[2293]298
[2299]299    CASE('tpot','q')
300      IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
301      ELSE;                 champ=champ*.01*workvar !--- Relative humidity
302        WHERE(champ<0.) champ=1.0E-10
303      END IF
304      DO il=1,lml
305        xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
306        xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
307        champ(:,1  ,il) = xppn
308        champ(:,jml,il) = xpps
309      END DO
310  END SELECT
[2293]311
312END SUBROUTINE startget_dyn3d
313!
314!-------------------------------------------------------------------------------
315
316
317!-------------------------------------------------------------------------------
318!
319SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
320!
321!-------------------------------------------------------------------------------
322  USE conf_dat_m, ONLY: conf_dat2d
323  IMPLICIT NONE
324!===============================================================================
325! Purpose:  Compute "phis" just like it would be in start_init_orog.
326!===============================================================================
327! Arguments:
328  REAL, INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
329  REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
330!-------------------------------------------------------------------------------
331! Local variables:
332  CHARACTER(LEN=256) :: modname="start_init_orog0"
333  CHARACTER(LEN=256) :: title="RELIEF"
334  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
335  REAL               :: lev(1), date, dt
336  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
337  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
338!-------------------------------------------------------------------------------
339  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
340  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
341  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
342  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
343  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
344  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
345
346!--- HIGH RESOLUTION OROGRAPHY
347  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
348
349  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
350  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
351                lev, ttm_tmp, itau, date, dt, fid)
352  ALLOCATE(relief_hi(iml_rel,jml_rel))
353  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
354  CALL flinclo(fid)
355
356!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
357  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
358  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
359  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
360
361!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
362  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
363  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
364  DEALLOCATE(lon_ini,lat_ini)
365
366!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
367  WRITE(lunout,*)
368  WRITE(lunout,*)'*** Compute surface geopotential ***'
369
370!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
371  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
372  phis = phis * 9.81
373  phis(iml,:) = phis(1,:)
374  DEALLOCATE(relief_hi,lon_rad,lat_rad)
375
376END SUBROUTINE start_init_orog0
377!
378!-------------------------------------------------------------------------------
379
380
381!-------------------------------------------------------------------------------
382!
383SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
384!
385!===============================================================================
386! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
387!          without any call to physics subroutines.
388!===============================================================================
389  IMPLICIT NONE
390!-------------------------------------------------------------------------------
391! Arguments:
392  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
393  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
394  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
395  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
396  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
397!-------------------------------------------------------------------------------
398! Local variables:
399  CHARACTER(LEN=256) :: modname="grid_noro0"
400  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
401  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
402  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
403  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
404  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
405  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
406  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
407  LOGICAL :: masque_lu
408  INTEGER :: i, ii, imdp, imar, iext
409  INTEGER :: j, jj, jmdp, jmar, nn
410  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
411  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
412!-------------------------------------------------------------------------------
413  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
414  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
415  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
416  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
417  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
418  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
419  iext=imdp/10
420  xpi = ACOS(-1.)
421  rad = 6371229.
422
423!--- ARE WE USING A READ MASK ?
424  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
425  WRITE(lunout,*)'Masque lu: ',masque_lu
426
427!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
428  ALLOCATE(xusn(imdp+2*iext))
429  xusn(1     +iext:imdp  +iext)=xd(:)
430  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
431  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
432
433  ALLOCATE(yusn(jmdp+2))
434  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
435  yusn(2:jmdp+1)=yd(:)
436  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
437
438  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
439  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
440  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
441  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
442  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
443  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
444  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
445  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
446
447!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
448  ALLOCATE(a(imar+1),b(imar+1))
449  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
450  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
451  a(1)=x(1)-(x(2)-x(1))/2.0
452  a(2:imar+1)= b(1:imar)
453
454  ALLOCATE(c(jmar),d(jmar))
455  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
456  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
457  c(1)=y(1)-(y(2)-y(1))/2.0
458  c(2:jmar)=d(1:jmar-1)
459
460!--- INITIALIZATIONS:
461  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
462  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
463
464!--- SUMMATION OVER GRIDPOINT AREA
465  zleny=xpi/REAL(jmdp)*rad
466  xincr=xpi/REAL(jmdp)/2.
467  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
468  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
469  DO ii = 1, imar+1
470    DO jj = 1, jmar
471      DO j = 2,jmdp+1
472        zlenx  =zleny  *COS(yusn(j))
473        zbordnor=(xincr+c(jj)-yusn(j))*rad
474        zbordsud=(xincr-d(jj)+yusn(j))*rad
475        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
476        IF(weighy/=0) THEN
477          DO i = 2, imdp+2*iext-1
478            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
479            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
480            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
481            IF(weighx/=0)THEN
482              num_tot(ii,jj)=num_tot(ii,jj)+1.0
483              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
484              weight(ii,jj)=weight(ii,jj)+weighx*weighy
485              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
486            END IF
487          END DO
488        END IF
489      END DO
490    END DO
491  END DO
492
493!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
494  IF(.NOT.masque_lu) THEN
495    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
496  END IF
497  nn=COUNT(weight(:,1:jmar-1)==0.0)
498  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
499  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
500
501!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
502  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
503  WHERE(mask>=0.1) mask_tmp = 1.
504  WHERE(weight(:,:)/=0.0)
505    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
506    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
507  END WHERE
508  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
509
510!--- Values at poles
511  zphi(imar+1,:)=zphi(1,:)
512
513  zweinor=SUM(weight(1:imar,   1),DIM=1)
514  zweisud=SUM(weight(1:imar,jmar),DIM=1)
515  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
516  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
517  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
518
519END SUBROUTINE grid_noro0
520!
521!-------------------------------------------------------------------------------
522
523
524!-------------------------------------------------------------------------------
525!
526SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol)
527!
528!-------------------------------------------------------------------------------
529  IMPLICIT NONE
530!===============================================================================
531! Purpose:   Compute psol, knowing phis.
532!===============================================================================
533! Arguments:
534  REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
535  REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
536  LOGICAL, INTENT(IN)  :: ibar
537  REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
538  REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
539!-------------------------------------------------------------------------------
540! Local variables:
541  CHARACTER(LEN=256) :: modname='start_init_dyn'
542  REAL               :: date, dt
543  INTEGER            :: iml, jml, jml2, itau(1)
544  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
545  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
546  REAL, ALLOCATABLE  :: z(:,:), ps(:,:), ts(:,:)
547!-------------------------------------------------------------------------------
548  iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2),            &
549      &                                              TRIM(modname)//" iml")
550  jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
551  jml2=SIZE(lat_in2)
552
553  WRITE(lunout,*) 'Opening the surface analysis'
554  CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
555  WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
556
557  ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
558  ALLOCATE(levdyn_ini(llm_dyn))
559  CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,                  &
560                lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn)
561
562!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
563  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
564  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
565  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
566
567  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
568  CALL get_var_dyn('Z',z)                        !--- SURFACE GEOPOTENTIAL
569  CALL get_var_dyn('SP',ps)                      !--- SURFACE PRESSURE
570  CALL get_var_dyn('ST',ts)                      !--- SURFACE TEMPERATURE
571!  CALL flinclo(fid_dyn)
572  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
573
574!--- PSOL IS COMPUTED IN PASCALS
575  psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0          &
576     &            /ts(:iml-1,:))
577  psol(iml,:)=psol(1,:)
578  DEALLOCATE(z,ps,ts)
579  psol(:,1  )=SUM(aire(1:iml-1,1  )*psol(1:iml-1,1  ))/apoln  !--- NORTH POLE
580  psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols  !--- SOUTH POLE
581
582CONTAINS
583
584!-------------------------------------------------------------------------------
585!
586SUBROUTINE get_var_dyn(title,field)
587!
588!-------------------------------------------------------------------------------
589  USE conf_dat_m, ONLY: conf_dat2d
590  IMPLICIT NONE
591!-------------------------------------------------------------------------------
592! Arguments:
593  CHARACTER(LEN=*),  INTENT(IN)    :: title
594  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
595!-------------------------------------------------------------------------------
596! Local variables:
597  CHARACTER(LEN=256) :: msg
598  INTEGER :: tllm
599!-------------------------------------------------------------------------------
600  SELECT CASE(title)
601    CASE('Z');     tllm=0;       msg='geopotential'
602    CASE('SP');    tllm=0;       msg='surface pressure'
603    CASE('ST');    tllm=llm_dyn; msg='temperature'
604  END SELECT
605  IF(.NOT.ALLOCATED(field)) THEN
606    ALLOCATE(field(iml,jml))
607    CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
608    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)
609    CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,       &
610                              lon_in, lat_in, lon_in2, lat_in2, field)
611  ELSE IF(SIZE(field)/=SIZE(z)) THEN
612    msg='The '//TRIM(msg)//' field we have does not have the right size'
613    CALL abort_gcm(TRIM(modname),msg,1)
614  END IF
615
616END SUBROUTINE get_var_dyn
617!
618!-------------------------------------------------------------------------------
619
620END SUBROUTINE start_init_dyn
621!
622!-------------------------------------------------------------------------------
623
624
625!-------------------------------------------------------------------------------
626!
627SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar)
628!
629!-------------------------------------------------------------------------------
630  USE conf_dat_m, ONLY: conf_dat3d
631  USE pchsp_95_m, ONLY: pchsp_95
632  USE pchfe_95_m, ONLY: pchfe_95
633  IMPLICIT NONE
634!-------------------------------------------------------------------------------
635! Arguments:
636  CHARACTER(LEN=*), INTENT(IN) :: var
637  REAL,    INTENT(IN)  :: lon_in(:),  lat_in(:)   ! dim (iml) (jml)
638  REAL,    INTENT(IN)  :: lon_in2(:), lat_in2(:)  ! dim (iml) (jml2)
639  REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
640  REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
641  LOGICAL, INTENT(IN)  :: ibar
642!-------------------------------------------------------------------------------
643! Local variables:
644  CHARACTER(LEN=256) :: modname='start_inter_3d'
645  LOGICAL :: skip
646  REAL    :: chmin, chmax
647  INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
648  INTEGER :: n_extrap                             ! Extrapolated points number
649  REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
650  REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
651  REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
652!-------------------------------------------------------------------------------
653  iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
654  jml=assert_eq(SIZE(lat_in),              SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
655  lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
656
657  WRITE(lunout, *)'Going into flinget to extract the 3D field.'
658  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
659  CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
660
661!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS
662  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
663  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
664  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
665
666!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
667  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
668  CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
669                       lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
670  DEALLOCATE(lon_ini, lat_ini)
671
672!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
673  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
674  DO il = 1,llm_dyn
675    CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),     &
676                          lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
677  END DO
678  DEALLOCATE(lon_rad, lat_rad)
679
680!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
681  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
682  ax = lev_dyn(llm_dyn:1:-1)
683  skip = .FALSE.
684  n_extrap = 0
685  DO ij=1, jml
686    DO ii=1, iml-1
687      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
688      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
689      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1),              &
690           var3d(ii, ij, lml:1:-1), ierr)
691      IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
692      n_extrap = n_extrap + ierr
693    END DO
694  END DO
695  IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
696  var3d(iml, :, :) = var3d(1, :, :)
697
698  DO il=1, lml
699    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
700    WRITE(lunout, *)' '//TRIM(var)//'  min max l ', il, chmin, chmax
701  END DO
702
703END SUBROUTINE start_inter_3d
704!
705!-------------------------------------------------------------------------------
706
707
708!-------------------------------------------------------------------------------
709!
710SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
711!
712!-------------------------------------------------------------------------------
713  USE inter_barxy_m, ONLY: inter_barxy
714  USE grid_atob_m,   ONLY: grille_m
715  IMPLICIT NONE
716!-------------------------------------------------------------------------------
717! Arguments:
718  CHARACTER(LEN=*), INTENT(IN)  :: nam
719  LOGICAL,          INTENT(IN)  :: ibar, ibeg
720  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
721  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
722  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
723  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
724  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
725!-------------------------------------------------------------------------------
726! Local variables:
727  CHARACTER(LEN=256) :: modname="interp_startvar"
728  INTEGER            :: ii, jj, i1, j1, j2
729  REAL, ALLOCATABLE  :: vtmp(:,:)
730!-------------------------------------------------------------------------------
731  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
732  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
733  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
734  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
735  j2=SIZE(lat2)
736  ALLOCATE(vtmp(i1-1,j1))
737  IF(ibar) THEN
738    IF(ibeg.AND.prt_level>1) THEN
739      WRITE(lunout,*)"---------------------------------------------------------"
740      WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
741      WRITE(lunout,*)"---------------------------------------------------------"
742    END IF
743    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
744  ELSE
745    CALL grille_m   (lon, lat,        vari, lon1,        lat1, vtmp)
746  END IF
747  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
748
749END SUBROUTINE interp_startvar
750!
751!-------------------------------------------------------------------------------
752
753END MODULE etat0dyn
754!
755!*******************************************************************************
Note: See TracBrowser for help on using the repository browser.