source: LMDZ6/branches/ICOLMDZISO/libf/phylmd/create_etat0_unstruct_mod.F90 @ 5594

Last change on this file since 5594 was 5592, checked in by yann meurdesoif, 8 months ago

Update ICOLMDZISO branch.
YM
File size: 11.0 KB
Line 
1MODULE create_etat0_unstruct_mod
2
3  REAL, SAVE, ALLOCATABLE :: zmea_gw(:)
4  !$OMP THREADPRIVATE(zmea_gw)
5  REAL, SAVE, ALLOCATABLE :: zpic_gw(:)
6  !$OMP THREADPRIVATE(zpic_gw)
7  REAL, SAVE, ALLOCATABLE :: zval_gw(:)
8  !$OMP THREADPRIVATE(zval_gw)
9  REAL, SAVE, ALLOCATABLE :: zstd_gw(:)
10  !$OMP THREADPRIVATE(zstd_gw)
11  REAL, SAVE, ALLOCATABLE :: zsig_gw(:)
12  !$OMP THREADPRIVATE(zsig_gw)
13  REAL, SAVE, ALLOCATABLE :: zgam_gw(:)
14  !$OMP THREADPRIVATE(zgam_gw)
15  REAL, SAVE, ALLOCATABLE :: zthe_gw(:)
16  !$OMP THREADPRIVATE(zthe_gw)
17
18  PRIVATE zmea_gw, zpic_gw, zval_gw, zstd_gw, zsig_gw, zgam_gw, zthe_gw
19
20
21CONTAINS
22 
23  SUBROUTINE init_create_etat0_unstruct
24  USE lmdz_xios
25  USE netcdf
26  USE mod_phys_lmdz_para
27  IMPLICIT NONE
28  INTEGER :: file_id, iret
29 
30   ! for coupling activate ocean fraction reading from file "ocean_fraction.nc"
31    IF (is_omp_master) THEN
32
33      IF (NF90_OPEN("ocean_fraction.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
34        CALL xios_set_file_attr("frac_ocean",enabled=.TRUE.)
35        CALL xios_set_field_attr("mask",field_ref="frac_ocean_read")
36        iret=NF90_CLOSE(file_id)
37      ELSE IF (NF90_OPEN("land_water_0.05.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
38        CALL xios_set_file_attr("land_water",name="land_water_0.05",enabled=.TRUE.)
39        CALL xios_set_field_attr("mask",field_ref="land_water")
40        iret=NF90_CLOSE(file_id)
41      ELSE IF (NF90_OPEN("land_water_0.25.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
42        CALL xios_set_file_attr("land_water",name="land_water_0.25",enabled=.TRUE.)
43        CALL xios_set_field_attr("mask",field_ref="land_water")
44        iret=NF90_CLOSE(file_id)
45      ELSE IF (NF90_OPEN("land_water_0.50.nc", NF90_NOWRITE, file_id)==NF90_NOERR) THEN
46        CALL xios_set_file_attr("land_water",name="land_water_0.50",enabled=.TRUE.)
47        CALL xios_set_field_attr("mask",field_ref="land_water")
48        iret=NF90_CLOSE(file_id)
49      ENDIF
50
51    ENDIF
52
53  END SUBROUTINE init_create_etat0_unstruct
54 
55 
56  SUBROUTINE init_param_gw(zmea, zpic, zval, zstd, zsig, zgam, zthe)
57  USE dimphy
58    REAL, INTENT(IN) :: zmea(klon) 
59    REAL, INTENT(IN) :: zpic(klon)
60    REAL, INTENT(IN) :: zval(klon)
61    REAL, INTENT(IN) :: zstd(klon)
62    REAL, INTENT(IN) :: zsig(klon)
63    REAL, INTENT(IN) :: zgam(klon)
64    REAL, INTENT(IN) :: zthe(klon)
65
66    ALLOCATE(zmea_gw(klon), zpic_gw(klon), zval_gw(klon), zstd_gw(klon), zsig_gw(klon), zgam_gw(klon), zthe_gw(klon))
67   
68    zmea_gw(:)=zmea(:)
69    zpic_gw(:)=zpic(:)
70    zval_gw(:)=zval(:)
71    zstd_gw(:)=zstd(:)
72    zsig_gw(:)=zsig(:)
73    zgam_gw(:)=zgam(:)
74    zthe_gw(:)=zthe(:)
75
76  END SUBROUTINE init_param_gw
77
78
79
80
81  SUBROUTINE create_etat0_unstruct
82  USE dimphy
83  USE lmdz_xios
84  USE infotrac_phy
85  USE fonte_neige_mod
86  USE pbl_surface_mod
87  USE phys_state_var_mod
88  USE indice_sol_mod
89  USE surface_data,      ONLY: landice_opt
90  USE mod_phys_lmdz_para
91  USE print_control_mod, ONLY: lunout
92  USE geometry_mod
93  USE ioipsl_getin_p_mod, ONLY: getin_p
94 
95  USE infotrac_phy, ONLY: niso
96  USE isotopes_routines_mod, ONLY: phyisoetat0
97  USE isotopes_mod, ONLY: iso_eau
98  USE isotopes_verif_mod, ONLY: iso_verif_egalite_vect2D,iso_verif_egalite
99
100  IMPLICIT NONE
101  INCLUDE 'dimsoil.h'
102  include "clesphys.h"
103
104    LOGICAL :: no_ter_antartique   ! If true, no land points are allowed at Antartic
105    REAL,    DIMENSION(klon)                 :: tsol
106    REAL,    DIMENSION(klon)                 :: sn
107    REAL,    DIMENSION(klon)                 :: rugmer
108    REAL,    DIMENSION(klon)                 :: run_off_lic_0
109    REAL,    DIMENSION(klon)                 :: lic
110    REAL,    DIMENSION(klon)                 :: fder
111
112    REAL,    DIMENSION(klon,nbsrf)           :: qsurf, snsrf
113    REAL,    DIMENSION(klon,nsoilmx,nbsrf)   :: tsoil
114   
115    REAL,    DIMENSION(klon_mpi)             :: tsol_mpi, qsol_mpi, zmasq_mpi, lic_mpi
116    REAL,    DIMENSION(klon_mpi)             :: zmea_mpi, zstd_mpi, zsig_mpi, zgam_mpi, zthe_mpi
117    REAL,    DIMENSION(klon_mpi)             :: cell_area_mpi
118    REAL,    DIMENSION(klon_mpi,nbsrf)       :: pctsrf_mpi
119#ifdef ISO
120    REAL xtsnow(niso,klon, nbsrf)
121    REAL xtrun_off_lic_0(niso,klon)
122    REAL Rland_ice(niso,klon)
123#endif
124
125    INCLUDE "compbl.h"
126    INCLUDE "alpale.h"
127   
128    INTEGER :: ji,j,i
129 
130
131!--- Initial atmospheric CO2 conc. from .def file
132    co2_ppm0 = co2_ppm
133
134    IF (is_omp_master) THEN
135      CALL xios_recv_field("ts",tsol_mpi)
136      CALL xios_recv_field("qs",qsol_mpi)
137      CALL xios_recv_field("mask",zmasq_mpi)
138      IF (landice_opt .LT. 2) CALL xios_recv_field("landice",lic_mpi)
139    ENDIF
140    CALL scatter_omp(tsol_mpi,tsol)
141    CALL scatter_omp(qsol_mpi,qsol)
142    CALL scatter_omp(zmasq_mpi,zmasq)
143    IF (landice_opt .LT. 2) CALL scatter_omp(lic_mpi,lic)
144
145    radsol(:)   = 0.0
146    rugmer(:) = 0.001
147    sn(:)     = 0
148
149    WHERE(qsol(:)<0) qsol(:)=0
150       
151    WHERE(   zmasq(:)<EPSFRA) zmasq(:)=0.
152    WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1.
153
154    pctsrf(:,:) = 0
155    IF (landice_opt .LT. 2) THEN
156       pctsrf(:,is_lic)=lic
157       WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
158       WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
159
160       pctsrf(:,is_ter)=zmasq(:)
161       
162       !--- Adequation with soil/sea mask
163       DO ji=1,klon
164          IF(zmasq(ji)>EPSFRA) THEN
165             IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
166                pctsrf(ji,is_lic)=zmasq(ji)
167                pctsrf(ji,is_ter)=0.
168             ELSE
169                pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
170                IF(pctsrf(ji,is_ter)<EPSFRA) THEN
171                   pctsrf(ji,is_ter)=0.
172                   pctsrf(ji,is_lic)=zmasq(ji)
173                END IF
174             END IF
175          END IF
176       END DO
177   
178    ELSE
179       ! landice_opt=>2 : no land ice
180       pctsrf(:,is_lic)=0.0
181       pctsrf(:,is_ter)=zmasq(:)
182    END IF
183
184
185
186
187
188  !--- Option no_ter_antartique removes all land fractions souther than 60S.
189  !--- Land ice is set instead of the land fractions on these latitudes.
190  !--- The ocean and sea-ice fractions are not changed.
191  !--- This option is only available if landice_opt<2.   
192  IF (landice_opt .LT. 2) THEN
193     no_ter_antartique=.FALSE.
194     CALL getin_p('no_ter_antartique',no_ter_antartique)
195     WRITE(lunout,*)"no_ter_antartique=",no_ter_antartique
196     IF (no_ter_antartique) THEN
197        ! Remove all land fractions souther than 60S and set land-ice instead
198        WRITE(lunout,*) "Remove land fractions souther than 60deg south by increasing"
199        WRITE(lunout,*) "the continental ice fractions. No land can now be found at Antartic."
200        DO ji=1, klon
201           IF (latitude_deg(ji)<-60.0) THEN
202              pctsrf(ji,is_lic) = pctsrf(ji,is_lic) + pctsrf(ji,is_ter)
203              pctsrf(ji,is_ter) = 0
204           END IF
205        END DO
206     END IF
207  END IF
208   
209! sub-surface ocean and sea ice (sea ice set to zero for start)
210!*******************************************************************************
211    pctsrf(:,is_oce)=(1.-zmasq(:))
212    WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
213   
214!    zval(:)=max(0.,zmea-2*zstd(:))
215!    zpic(:)=zmea+2*zstd(:)
216   
217!! WARNING    DON'T FORGET FOR LATER
218!!ym  IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
219!!
220   
221! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
222!*******************************************************************************
223    DO i=1,nbsrf
224     ftsol(:,i) = tsol
225    END DO
226 
227    DO i=1,nbsrf
228     snsrf(:,i) = sn
229    END DO
230!albedo SB >>>
231!ym error : the sub surface dimension is the third not second
232!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
233!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
234    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
235    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
236
237!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
238!ym probably the uninitialized value was 0 for standard (regular grid) case
239    falb_dif(:,:,:)=0
240    u10m(:,:)=0 
241    v10m(:,:)=0 
242    treedrg(:,:,:)=0
243!albedo SB <<<
244    fevap(:,:) = 0.
245    qsurf = 0.
246   
247    DO i=1,nbsrf
248      DO j=1,nsoilmx
249        tsoil(:,j,i) = tsol
250      END DO
251    END DO
252 
253    rain_fall = 0.; snow_fall = 0.
254    solsw = 165.;   sollw = -53.
255    solswfdiff = 1.
256!ym warning missing init for sollwdown => set to 0
257  sollwdown  = 0.
258   
259   
260    t_ancien = 273.15
261    u_ancien=0
262    v_ancien=0
263    q_ancien = 0.
264    ql_ancien = 0.
265    qs_ancien = 0.
266    prlw_ancien = 0.
267    prsw_ancien = 0.
268    prw_ancien = 0.
269    agesno = 0.
270   
271    ! LSCP condensation and ice supersaturation
272    cf_ancien = 0.
273    rvc_ancien = 0.
274
275    wake_delta_pbl_TKE(:,:,:)=0
276    wake_dens(:)=0
277    awake_dens = 0.
278    cv_gen = 0.
279    ale_bl = 0.
280    ale_bl_trig =0.
281    alp_bl=0.
282    ale_wake=0.
283    ale_bl_stat=0.
284   
285    z0m(:,:)=0 ! ym missing 5th subsurface initialization
286    z0m(:,is_oce) = rugmer(:)
287    z0m(:,is_ter) = 0.01 ! MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
288    z0m(:,is_lic) = 0.001 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
289    z0m(:,is_sic) = 0.001
290    z0h(:,:)=z0m(:,:)
291
292    fder = 0.0
293    clwcon = 0.0
294    rnebcon = 0.0
295    ratqs = 0.0
296    run_off_lic_0 = 0.0
297    rugoro = 0.0
298
299! Before phyredem calling, surface modules and values to be saved in startphy.nc
300! are initialized
301!*******************************************************************************
302    pbl_tke(:,:,:) = 1.e-8
303    zmax0(:) = 40.
304    f0(:) = 1.e-5
305    sig1(:,:) = 0.
306    w01(:,:) = 0.
307    wake_deltat(:,:) = 0.
308    wake_deltaq(:,:) = 0.
309    wake_s(:) = 0.
310    wake_cstar(:) = 0.
311    wake_fip(:) = 0.
312    wake_pe = 0.
313    fm_therm = 0.
314    entr_therm = 0.
315    detr_therm = 0.
316    awake_s = 0.
317
318
319    CALL fonte_neige_init(run_off_lic_0)
320    CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
321
322    IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then
323     delta_tsurf = 0.
324     beta_aridity = 0.
325    END IF
326    ratqs_inter_ = 0.002
327
328
329#ifdef ISO
330        ! initialise les isotopes       
331        write(*,*) 'phyetat0 1069'
332         CALL phyisoetat0 (snow,run_off_lic_0, &
333     &           xtsnow,xtrun_off_lic_0, &
334     &           Rland_ice)
335#ifdef ISOVERIF
336      write(*,*) 'phyetat0 1074'
337      if (iso_eau.gt.0) then
338      call iso_verif_egalite_vect2D(  &
339     &           xtsnow,snow, &
340     &           'phyetat0 1101a',niso,klon,nbsrf)
341        do i=1,klon 
342              call iso_verif_egalite(Rland_ice(iso_eau,i),1.0, &
343     &         'phyetat0 1101b')
344         enddo
345      endif
346      write(*,*) 'phyetat0 1102'
347#endif
348#endif
349
350#ifdef ISO
351   CALL fonte_neige_init_iso(xtrun_off_lic_0)
352#endif
353
354
355#ifdef ISO
356  CALL pbl_surface_init_iso(xtsnow,Rland_ice)
357#endif
358
359    CALL gather_omp(cell_area,cell_area_mpi)
360    CALL gather_omp(pctsrf,pctsrf_mpi)
361    IF (is_omp_master) THEN
362      CALL xios_send_field("area_ce0l",cell_area_mpi)
363      CALL xios_send_field("fract_oce_ce0l",pctsrf_mpi(:,is_oce))
364      CALL xios_send_field("fract_sic_ce0l",pctsrf_mpi(:,is_sic))
365    ENDIF
366   
367    zmea(:) = zmea_gw(:)
368    zpic(:) = zpic_gw(:)
369    zval(:) = zval_gw(:)
370    zstd(:) = zstd_gw(:)
371    zsig(:) = zsig_gw(:)
372    zgam(:) = zgam_gw(:)
373    zthe(:) = zthe_gw(:)
374    DEALLOCATE(zmea_gw, zpic_gw, zval_gw, zstd_gw, zsig_gw, zgam_gw, zthe_gw)
375
376    CALL phyredem( "startphy.nc" )
377
378  END SUBROUTINE create_etat0_unstruct
379
380
381END MODULE create_etat0_unstruct_mod
Note: See TracBrowser for help on using the repository browser.