source: LMDZ5/branches/AI-cosp/libf/phylmd/phyaqua_mod.F90 @ 5444

Last change on this file since 5444 was 2351, checked in by Ehouarn Millour, 9 years ago

More on physics/dynamics separation and cleanup:

  • Set things up so that all physics-related initializations are done via iniphysiq.
  • Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
  • removed obsolete comgeomphy.h and comgeomphy.F90

EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 28.9 KB
Line 
1MODULE phyaqua_mod
2  ! Routines complementaires pour la physique planetaire.
3  IMPLICIT NONE
4
5CONTAINS
6
7  SUBROUTINE iniaqua(nlon, iflag_phys)
8
9    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10    ! Creation d'un etat initial et de conditions aux limites
11    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
12    ! du modele LMDZ dans sa version terrestre.
13    ! iflag_phys est un parametre qui controle
14    ! iflag_phys = N
15    ! de 100 a 199 : aqua planetes avec SST forcees
16    ! N-100 determine le type de SSTs
17    ! de 200 a 299 : terra planetes avec Ts calcule
18
19    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20
21    USE dimphy, ONLY: klon
22    USE geometry_mod, ONLY : latitude
23    USE surface_data, ONLY: type_ocean, ok_veget
24    USE pbl_surface_mod, ONLY: pbl_surface_init
25    USE fonte_neige_mod, ONLY: fonte_neige_init
26    USE phys_state_var_mod
27    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
28                                day_ini,day_end
29    USE indice_sol_mod
30    USE nrtype, ONLY: pi
31    USE ioipsl
32    IMPLICIT NONE
33
34    include "YOMCST.h"
35    include "clesphys.h"
36    include "dimsoil.h"
37
38    INTEGER, INTENT (IN) :: nlon, iflag_phys
39    ! IM ajout latfi, lonfi
40!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
41
42    INTEGER type_profil, type_aqua
43
44    ! Ajouts initialisation des surfaces
45    REAL :: run_off_lic_0(nlon)
46    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
47    REAL :: tsoil(nlon, nsoilmx, nbsrf)
48    REAL :: tslab(nlon), seaice(nlon)
49    REAL fder(nlon)
50
51
52
53    ! Arguments :
54    ! -----------
55
56    ! integer radpas
57    INTEGER it, unit, i, k, itap
58
59    REAL airefi, zcufi, zcvfi
60
61    REAL rugos, albedo
62    REAL tsurf
63    REAL time, timestep, day, day0
64    REAL qsol_f
65    REAL rugsrel(nlon)
66    ! real zmea(nlon),zstd(nlon),zsig(nlon)
67    ! real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
68    ! real rlon(nlon),rlat(nlon)
69    LOGICAL alb_ocean
70    ! integer demih_pas
71
72    CHARACTER *80 ans, file_forctl, file_fordat, file_start
73    CHARACTER *100 file, var
74    CHARACTER *2 cnbl
75
76    REAL phy_nat(nlon, 360)
77    REAL phy_alb(nlon, 360)
78    REAL phy_sst(nlon, 360)
79    REAL phy_bil(nlon, 360)
80    REAL phy_rug(nlon, 360)
81    REAL phy_ice(nlon, 360)
82    REAL phy_fter(nlon, 360)
83    REAL phy_foce(nlon, 360)
84    REAL phy_fsic(nlon, 360)
85    REAL phy_flic(nlon, 360)
86
87    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
88
89    ! intermediate variables to use getin (need to be "save" to be shared by
90    ! all threads)
91    INTEGER, SAVE :: nbapp_rad_omp
92    REAL, SAVE :: co2_ppm_omp, solaire_omp
93    LOGICAL, SAVE :: alb_ocean_omp
94    REAL, SAVE :: rugos_omp
95    ! -------------------------------------------------------------------------
96    ! declaration pour l'appel a phyredem
97    ! -------------------------------------------------------------------------
98
99    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
100    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
101    ! real pbl_tke(nlon,llm,nbsrf)
102    ! real rain_fall(nlon),snow_fall(nlon)
103    ! real solsw(nlon), sollw(nlon),radsol(nlon)
104    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
105    ! real ratqs(nlon,llm)
106    ! real clwcon(nlon,llm)
107
108    INTEGER longcles
109    PARAMETER (longcles=20)
110    REAL clesphy0(longcles)
111
112
113    ! -----------------------------------------------------------------------
114    ! dynamial tendencies :
115    ! ---------------------
116
117    INTEGER l, ierr, aslun
118
119!    REAL longitude, latitude
120    REAL paire
121
122!    DATA latitude, longitude/48., 0./
123
124    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
125    ! INITIALISATIONS
126    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
127
128    ! -----------------------------------------------------------------------
129    ! Initialisations  des constantes
130    ! -------------------------------
131
132
133    type_aqua = iflag_phys/100
134    type_profil = iflag_phys - type_aqua*100
135    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
136
137    IF (klon/=nlon) THEN
138      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
139      STOP 'probleme de dimensions dans iniaqua'
140    END IF
141    CALL phys_state_var_init(read_climoz)
142
143
144    read_climoz = 0
145    day0 = 217.
146    day = day0
147    it = 0
148    time = 0.
149
150    ! IM ajout latfi, lonfi
151!    rlatd = latfi
152!    rlond = lonfi
153!    rlat = rlatd*180./pi
154!    rlon = rlond*180./pi
155
156    ! -----------------------------------------------------------------------
157    ! initialisations de la physique
158    ! -----------------------------------------------------------------------
159
160    day_ini = day_ref
161    day_end = day_ini + ndays
162!    airefi = 1.
163!    zcufi = 1.
164!    zcvfi = 1.
165    !$OMP MASTER
166    nbapp_rad_omp = 24
167    CALL getin('nbapp_rad', nbapp_rad_omp)
168    !$OMP END MASTER
169    !$OMP BARRIER
170    nbapp_rad = nbapp_rad_omp
171
172    ! ---------------------------------------------------------------------
173    ! Creation des conditions aux limites:
174    ! ------------------------------------
175    ! Initialisations des constantes
176    ! Ajouter les manquants dans planete.def... (albedo etc)
177    !$OMP MASTER
178    co2_ppm_omp = 348.
179    CALL getin('co2_ppm', co2_ppm_omp)
180    solaire_omp = 1365.
181    CALL getin('solaire', solaire_omp)
182    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
183    ! type_aqua
184    alb_ocean_omp = .TRUE.
185    CALL getin('alb_ocean', alb_ocean_omp)
186    !$OMP END MASTER
187    !$OMP BARRIER
188    co2_ppm = co2_ppm_omp
189    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
190    solaire = solaire_omp
191    WRITE (*, *) 'iniaqua: solaire=', solaire
192    alb_ocean = alb_ocean_omp
193    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
194
195    radsol = 0.
196    qsol_f = 10.
197
198    ! Conditions aux limites:
199    ! -----------------------
200
201    qsol(:) = qsol_f
202    rugsrel = 0.0 ! (rugsrel = rugoro)
203    rugoro = 0.0
204    u_ancien = 0.0
205    v_ancien = 0.0
206    agesno = 50.0
207    ! Relief plat
208    zmea = 0.
209    zstd = 0.
210    zsig = 0.
211    zgam = 0.
212    zthe = 0.
213    zpic = 0.
214    zval = 0.
215
216    ! Une seule surface
217    pctsrf = 0.
218    IF (type_aqua==1) THEN
219      rugos = 1.E-4
220      albedo = 0.19
221      pctsrf(:, is_oce) = 1.
222    ELSE IF (type_aqua==2) THEN
223      rugos = 0.03
224      albedo = 0.1
225      pctsrf(:, is_ter) = 1.
226    END IF
227
228    !$OMP MASTER
229    rugos_omp = rugos
230    CALL getin('rugos', rugos_omp)
231    !$OMP END MASTER
232    !$OMP BARRIER
233    rugos = rugos_omp
234    WRITE (*, *) 'iniaqua: rugos=', rugos
235    zmasq(:) = pctsrf(:, is_ter)
236
237    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
238    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
239
240    ! Si alb_ocean on calcule un albedo oceanique moyen
241    ! if (alb_ocean) then
242    ! Voir pourquoi on avait ca.
243    ! CALL ini_alb_oce(phy_alb)
244    ! else
245    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
246    ! endif !alb_ocean
247
248    DO i = 1, 360
249      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
250      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
251
252      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
253      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
254      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
255      phy_ice(:, i) = 0.0 ! fraction de glace (?)
256      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
257      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
258      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
259      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
260    END DO
261    ! IM calcul profil sst
262    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
263
264    CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
265      phy_fter, phy_foce, phy_flic, phy_fsic)
266
267
268    ! ---------------------------------------------------------------------
269    ! Ecriture de l'etat initial:
270    ! ---------------------------
271
272
273    ! Ecriture etat initial physique
274
275    timestep = pdtphys
276    radpas = nint(rday/timestep/float(nbapp_rad))
277
278    DO i = 1, longcles
279      clesphy0(i) = 0.
280    END DO
281    clesphy0(1) = float(iflag_con)
282    clesphy0(2) = float(nbapp_rad)
283    ! IF( cycle_diurne  ) clesphy0(3) =  1.
284    clesphy0(3) = 1. ! cycle_diurne
285    clesphy0(4) = 1. ! soil_model
286    clesphy0(5) = 1. ! new_oliq
287    clesphy0(6) = 0. ! ok_orodr
288    clesphy0(7) = 0. ! ok_orolf
289    clesphy0(8) = 0. ! ok_limitvrai
290
291
292    ! =======================================================================
293    ! Profils initiaux
294    ! =======================================================================
295
296    ! On initialise les temperatures de surfaces comme les sst
297    DO i = 1, nlon
298      ftsol(i, :) = phy_sst(i, 1)
299      tsoil(i, :, :) = phy_sst(i, 1)
300      tslab(i) = phy_sst(i, 1)
301    END DO
302
303    falbe(:, :) = albedo
304    falblw(:, :) = albedo
305    rain_fall(:) = 0.
306    snow_fall(:) = 0.
307    solsw(:) = 0.
308    sollw(:) = 0.
309    radsol(:) = 0.
310
311    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
312    ! intialisation bidon mais pas grave
313    t_ancien(:, :) = 0.
314    q_ancien(:, :) = 0.
315    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
316    rnebcon = 0.
317    ratqs = 0.
318    clwcon = 0.
319    pbl_tke = 1.E-8
320
321    ! variables supplementaires pour appel a plb_surface_init
322    fder(:) = 0.
323    seaice(:) = 0.
324    run_off_lic_0 = 0.
325    fevap = 0.
326
327
328    ! Initialisations necessaires avant phyredem
329    type_ocean = 'force'
330    CALL fonte_neige_init(run_off_lic_0)
331    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
332    snsrf(:, :) = 0. ! couverture de neige des sous surface
333    z0m(:, :) = rugos ! couverture de neige des sous surface
334    z0h=z0m
335
336
337    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
338
339    PRINT *, 'iniaqua: before phyredem'
340
341    falb1 = albedo
342    falb2 = albedo
343    zmax0 = 0.
344    f0 = 0.
345    sig1 = 0.
346    w01 = 0.
347    wake_deltat = 0.
348    wake_deltaq = 0.
349    wake_s = 0.
350    wake_cstar = 0.
351    wake_pe = 0.
352    wake_fip = 0.
353    fm_therm = 0.
354    entr_therm = 0.
355    detr_therm = 0.
356
357
358    CALL phyredem('startphy.nc')
359
360    PRINT *, 'iniaqua: after phyredem'
361    CALL phys_state_var_end
362
363    RETURN
364  END SUBROUTINE iniaqua
365
366
367  ! ====================================================================
368  ! ====================================================================
369  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
370    USE dimphy
371    IMPLICIT NONE
372    ! ====================================================================
373    ! =============================================================
374    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
375    ! Auteur : A. Campoy et F. Hourdin
376    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
377    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
378    ! connaissant la declinaison, la latitude et la longitude.
379
380    ! Dans cette version particuliere, on calcule le rayonnement
381    ! moyen sur l'année à chaque latitude.
382    ! angle zenithal calculé pour obtenir un
383    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
384    ! en moyenne annuelle.
385    ! Spécifique de la terre. Utilisé pour les aqua planetes.
386
387    ! Rque   : Different de la routine angle en ce sens que zenang
388    ! fournit des moyennes de pmu0 et non des valeurs
389    ! instantanees, du coup frac prend toutes les valeurs
390    ! entre 0 et 1.
391    ! Date   : premiere version le 13 decembre 1994
392    ! revu pour  GCM  le 30 septembre 1996
393    ! ===============================================================
394    ! longi----INPUT : la longitude vraie de la terre dans son plan
395    ! solaire a partir de l'equinoxe de printemps (degre)
396    ! gmtime---INPUT : temps universel en fraction de jour
397    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
398    ! lat------INPUT : latitude en degres
399    ! long-----INPUT : longitude en degres
400    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
401    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
402    ! ================================================================
403    include "YOMCST.h"
404    ! ================================================================
405    LOGICAL cycle_diurne
406    REAL gmtime
407    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
408    ! ================================================================
409    INTEGER i
410    REAL gmtime1, gmtime2
411    REAL pi_local
412
413
414    REAL rmu0m(klon), rmu0a(klon)
415
416
417    pi_local = 4.0*atan(1.0)
418
419    ! ================================================================
420    ! Calcul de l'angle zenithal moyen sur la journee
421    ! ================================================================
422
423    DO i = 1, klon
424      fract(i) = 1.
425      ! Calcule du flux moyen
426      IF (abs(rlat(i))<=28.75) THEN
427        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
428      ELSE IF (abs(rlat(i))<=43.75) THEN
429        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
430      ELSE IF (abs(rlat(i))<=71.25) THEN
431        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
432      ELSE
433        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
434      END IF
435    END DO
436
437    ! ================================================================
438    ! Avec ou sans cycle diurne
439    ! ================================================================
440
441    IF (cycle_diurne) THEN
442
443      ! On redecompose flux  au sommet suivant un cycle diurne idealise
444      ! identique a toutes les latitudes.
445
446      DO i = 1, klon
447        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
448        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
449          rmu0a(i)/sqrt(2.)
450      END DO
451
452      DO i = 1, klon
453        IF (rmu0(i)<=0.) THEN
454          rmu0(i) = 0.
455          fract(i) = 0.
456        ELSE
457          fract(i) = 1.
458        END IF
459      END DO
460
461      ! Affichage de l'angel zenitale
462      ! print*,'************************************'
463      ! print*,'************************************'
464      ! print*,'************************************'
465      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
466      ! print*,'rmu0m=',rmu0m(i)
467      ! print*,'rmu0a=',rmu0a(i)
468      ! print*,'rmu0=',rmu0(i)
469
470    ELSE
471
472      DO i = 1, klon
473        fract(i) = 0.5
474        rmu0(i) = rmu0m(i)*2.
475      END DO
476
477    END IF
478
479    RETURN
480  END SUBROUTINE zenang_an
481
482  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483
484  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
485      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
486
487    USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_omp_root
488    USE mod_grid_phy_lmdz, ONLY: klon_glo
489    USE mod_phys_lmdz_transfert_para, ONLY: gather
490    IMPLICIT NONE
491    include "netcdf.inc"
492
493    INTEGER, INTENT (IN) :: klon
494    REAL, INTENT (IN) :: phy_nat(klon, 360)
495    REAL, INTENT (IN) :: phy_alb(klon, 360)
496    REAL, INTENT (IN) :: phy_sst(klon, 360)
497    REAL, INTENT (IN) :: phy_bil(klon, 360)
498    REAL, INTENT (IN) :: phy_rug(klon, 360)
499    REAL, INTENT (IN) :: phy_ice(klon, 360)
500    REAL, INTENT (IN) :: phy_fter(klon, 360)
501    REAL, INTENT (IN) :: phy_foce(klon, 360)
502    REAL, INTENT (IN) :: phy_flic(klon, 360)
503    REAL, INTENT (IN) :: phy_fsic(klon, 360)
504
505    REAL :: phy_glo(klon_glo, 360) ! temporary variable, to store phy_***(:)
506      ! on the whole physics grid
507    INTEGER :: k
508    INTEGER ierr
509    INTEGER dimfirst(3)
510    INTEGER dimlast(3)
511
512    INTEGER nid, ndim, ntim
513    INTEGER dims(2), debut(2), epais(2)
514    INTEGER id_tim
515    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
516    INTEGER id_fter, id_foce, id_fsic, id_flic
517
518    IF (is_mpi_root .AND. is_omp_root) THEN
519
520      PRINT *, 'writelim: Ecriture du fichier limit'
521
522      ierr = nf_create('limit.nc', nf_clobber, nid)
523
524      ierr = nf_put_att_text(nid, nf_global, 'title', 30, &
525        'Fichier conditions aux limites')
526      ! !        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
527      ierr = nf_def_dim(nid, 'points_physiques', klon_glo, ndim)
528      ierr = nf_def_dim(nid, 'time', nf_unlimited, ntim)
529
530      dims(1) = ndim
531      dims(2) = ntim
532
533#ifdef NC_DOUBLE
534      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
535#else
536      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
537#endif
538      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
539
540#ifdef NC_DOUBLE
541      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
542#else
543      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
544#endif
545      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
546        'Nature du sol (0,1,2,3)')
547
548#ifdef NC_DOUBLE
549      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
550#else
551      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
552#endif
553      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
554        'Temperature superficielle de la mer')
555
556#ifdef NC_DOUBLE
557      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
558#else
559      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
560#endif
561      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
562        'Reference flux de chaleur au sol')
563
564#ifdef NC_DOUBLE
565      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
566#else
567      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
568#endif
569      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
570
571#ifdef NC_DOUBLE
572      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
573#else
574      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
575#endif
576      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
577
578#ifdef NC_DOUBLE
579      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
580#else
581      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
582#endif
583      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
584#ifdef NC_DOUBLE
585      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
586#else
587      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
588#endif
589      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
590#ifdef NC_DOUBLE
591      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
592#else
593      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
594#endif
595      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
596#ifdef NC_DOUBLE
597      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
598#else
599      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
600#endif
601      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
602
603      ierr = nf_enddef(nid)
604      IF (ierr/=nf_noerr) THEN
605        WRITE (*, *) 'writelim error: failed to end define mode'
606        WRITE (*, *) nf_strerror(ierr)
607      END IF
608
609
610      ! write the 'times'
611      DO k = 1, 360
612#ifdef NC_DOUBLE
613        ierr = nf_put_var1_double(nid, id_tim, k, dble(k))
614#else
615        ierr = nf_put_var1_real(nid, id_tim, k, float(k))
616#endif
617        IF (ierr/=nf_noerr) THEN
618          WRITE (*, *) 'writelim error with temps(k),k=', k
619          WRITE (*, *) nf_strerror(ierr)
620        END IF
621      END DO
622
623    END IF ! of if (is_mpi_root.and.is_omp_root)
624
625    ! write the fields, after having collected them on master
626
627    CALL gather(phy_nat, phy_glo)
628    IF (is_mpi_root .AND. is_omp_root) THEN
629#ifdef NC_DOUBLE
630      ierr = nf_put_var_double(nid, id_nat, phy_glo)
631#else
632      ierr = nf_put_var_real(nid, id_nat, phy_glo)
633#endif
634      IF (ierr/=nf_noerr) THEN
635        WRITE (*, *) 'writelim error with phy_nat'
636        WRITE (*, *) nf_strerror(ierr)
637      END IF
638    END IF
639
640    CALL gather(phy_sst, phy_glo)
641    IF (is_mpi_root .AND. is_omp_root) THEN
642#ifdef NC_DOUBLE
643      ierr = nf_put_var_double(nid, id_sst, phy_glo)
644#else
645      ierr = nf_put_var_real(nid, id_sst, phy_glo)
646#endif
647      IF (ierr/=nf_noerr) THEN
648        WRITE (*, *) 'writelim error with phy_sst'
649        WRITE (*, *) nf_strerror(ierr)
650      END IF
651    END IF
652
653    CALL gather(phy_bil, phy_glo)
654    IF (is_mpi_root .AND. is_omp_root) THEN
655#ifdef NC_DOUBLE
656      ierr = nf_put_var_double(nid, id_bils, phy_glo)
657#else
658      ierr = nf_put_var_real(nid, id_bils, phy_glo)
659#endif
660      IF (ierr/=nf_noerr) THEN
661        WRITE (*, *) 'writelim error with phy_bil'
662        WRITE (*, *) nf_strerror(ierr)
663      END IF
664    END IF
665
666    CALL gather(phy_alb, phy_glo)
667    IF (is_mpi_root .AND. is_omp_root) THEN
668#ifdef NC_DOUBLE
669      ierr = nf_put_var_double(nid, id_alb, phy_glo)
670#else
671      ierr = nf_put_var_real(nid, id_alb, phy_glo)
672#endif
673      IF (ierr/=nf_noerr) THEN
674        WRITE (*, *) 'writelim error with phy_alb'
675        WRITE (*, *) nf_strerror(ierr)
676      END IF
677    END IF
678
679    CALL gather(phy_rug, phy_glo)
680    IF (is_mpi_root .AND. is_omp_root) THEN
681#ifdef NC_DOUBLE
682      ierr = nf_put_var_double(nid, id_rug, phy_glo)
683#else
684      ierr = nf_put_var_real(nid, id_rug, phy_glo)
685#endif
686      IF (ierr/=nf_noerr) THEN
687        WRITE (*, *) 'writelim error with phy_rug'
688        WRITE (*, *) nf_strerror(ierr)
689      END IF
690    END IF
691
692    CALL gather(phy_fter, phy_glo)
693    IF (is_mpi_root .AND. is_omp_root) THEN
694#ifdef NC_DOUBLE
695      ierr = nf_put_var_double(nid, id_fter, phy_glo)
696#else
697      ierr = nf_put_var_real(nid, id_fter, phy_glo)
698#endif
699      IF (ierr/=nf_noerr) THEN
700        WRITE (*, *) 'writelim error with phy_fter'
701        WRITE (*, *) nf_strerror(ierr)
702      END IF
703    END IF
704
705    CALL gather(phy_foce, phy_glo)
706    IF (is_mpi_root .AND. is_omp_root) THEN
707#ifdef NC_DOUBLE
708      ierr = nf_put_var_double(nid, id_foce, phy_glo)
709#else
710      ierr = nf_put_var_real(nid, id_foce, phy_glo)
711#endif
712      IF (ierr/=nf_noerr) THEN
713        WRITE (*, *) 'writelim error with phy_foce'
714        WRITE (*, *) nf_strerror(ierr)
715      END IF
716    END IF
717
718    CALL gather(phy_fsic, phy_glo)
719    IF (is_mpi_root .AND. is_omp_root) THEN
720#ifdef NC_DOUBLE
721      ierr = nf_put_var_double(nid, id_fsic, phy_glo)
722#else
723      ierr = nf_put_var_real(nid, id_fsic, phy_glo)
724#endif
725      IF (ierr/=nf_noerr) THEN
726        WRITE (*, *) 'writelim error with phy_fsic'
727        WRITE (*, *) nf_strerror(ierr)
728      END IF
729    END IF
730
731    CALL gather(phy_flic, phy_glo)
732    IF (is_mpi_root .AND. is_omp_root) THEN
733#ifdef NC_DOUBLE
734      ierr = nf_put_var_double(nid, id_flic, phy_glo)
735#else
736      ierr = nf_put_var_real(nid, id_flic, phy_glo)
737#endif
738      IF (ierr/=nf_noerr) THEN
739        WRITE (*, *) 'writelim error with phy_flic'
740        WRITE (*, *) nf_strerror(ierr)
741      END IF
742    END IF
743
744    ! close file:
745    IF (is_mpi_root .AND. is_omp_root) THEN
746      ierr = nf_close(nid)
747    END IF
748
749  END SUBROUTINE writelim
750
751  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
752
753  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
754    USE dimphy
755    IMPLICIT NONE
756
757    INTEGER nlon, type_profil, i, k, j
758    REAL :: rlatd(nlon), phy_sst(nlon, 360)
759    INTEGER imn, imx, amn, amx, kmn, kmx
760    INTEGER p, pplus, nlat_max
761    PARAMETER (nlat_max=72)
762    REAL x_anom_sst(nlat_max)
763
764    IF (klon/=nlon) STOP 'probleme de dimensions dans iniaqua'
765    WRITE (*, *) ' profil_sst: type_profil=', type_profil
766    DO i = 1, 360
767      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
768
769      ! Rajout fbrlmd
770
771      IF (type_profil==1) THEN
772        ! Méthode 1 "Control" faible plateau à l'Equateur
773        DO j = 1, klon
774          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
775          ! PI/3=1.047197551
776          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
777            phy_sst(j, i) = 273.
778          END IF
779        END DO
780      END IF
781      IF (type_profil==2) THEN
782        ! Méthode 2 "Flat" fort plateau à l'Equateur
783        DO j = 1, klon
784          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
785          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
786            phy_sst(j, i) = 273.
787          END IF
788        END DO
789      END IF
790
791
792      IF (type_profil==3) THEN
793        ! Méthode 3 "Qobs" plateau réel à l'Equateur
794        DO j = 1, klon
795          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
796            rlatd(j))**4)
797          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
798            phy_sst(j, i) = 273.
799          END IF
800        END DO
801      END IF
802
803      IF (type_profil==4) THEN
804        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
805        DO j = 1, klon
806          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
807            rlatd(j))**4)
808          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
809            phy_sst(j, i) = 273.
810          END IF
811        END DO
812      END IF
813
814      IF (type_profil==5) THEN
815        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
816        DO j = 1, klon
817          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
818            *rlatd(j))**4)
819          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
820            phy_sst(j, i) = 273. + 2.
821          END IF
822
823        END DO
824      END IF
825
826      IF (type_profil==6) THEN
827        ! Méthode 6 "cst" valeur constante de SST
828        DO j = 1, klon
829          phy_sst(j, i) = 288.
830        END DO
831      END IF
832
833
834      IF (type_profil==7) THEN
835        ! Méthode 7 "cst" valeur constante de SST +2
836        DO j = 1, klon
837          phy_sst(j, i) = 288. + 2.
838        END DO
839      END IF
840
841      p = 0
842      IF (type_profil==8) THEN
843        ! Méthode 8 profil anomalies SST du modèle couplé AR4
844        DO j = 1, klon
845          IF (rlatd(j)==rlatd(j-1)) THEN
846            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
847              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
848            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
849              phy_sst(j, i) = 273. + x_anom_sst(pplus)
850            END IF
851          ELSE
852            p = p + 1
853            pplus = 73 - p
854            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
855              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
856            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
857              phy_sst(j, i) = 273. + x_anom_sst(pplus)
858            END IF
859            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
860          END IF
861        END DO
862      END IF
863
864      IF (type_profil==9) THEN
865        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
866        DO j = 1, klon
867          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
868            *rlatd(j))**4)
869          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
870            phy_sst(j, i) = 273. - 2.
871          END IF
872        END DO
873      END IF
874
875
876      IF (type_profil==10) THEN
877        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
878        DO j = 1, klon
879          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
880            *rlatd(j))**4)
881          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
882            phy_sst(j, i) = 273. + 4.
883          END IF
884        END DO
885      END IF
886
887      IF (type_profil==11) THEN
888        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
889        DO j = 1, klon
890          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
891            rlatd(j))**4)
892          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
893            phy_sst(j, i) = 273.
894          END IF
895        END DO
896      END IF
897
898      IF (type_profil==12) THEN
899        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
900        DO j = 1, klon
901          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
902            *rlatd(j))**4)
903          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
904            phy_sst(j, i) = 273. + 4.
905          END IF
906        END DO
907      END IF
908
909      IF (type_profil==13) THEN
910        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
911        DO j = 1, klon
912          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
913            rlatd(j))**4)
914          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
915            phy_sst(j, i) = 273.
916          END IF
917        END DO
918      END IF
919
920      IF (type_profil==14) THEN
921        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
922        DO j = 1, klon
923          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
924            *rlatd(j))**4)
925          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
926            phy_sst(j, i) = 273.
927          END IF
928        END DO
929      END IF
930
931      if (type_profil.EQ.20) then
932      print*,'Profile SST 20'
933!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
934
935      do j=1,klon
936        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
937      enddo
938      endif
939
940      if (type_profil.EQ.21) then
941      print*,'Profile SST 21'
942!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
943      do j=1,klon
944        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
945      enddo
946      endif
947
948
949
950    END DO
951
952    ! IM beg : verif profil SST: phy_sst
953    amn = min(phy_sst(1,1), 1000.)
954    amx = max(phy_sst(1,1), -1000.)
955    imn = 1
956    kmn = 1
957    imx = 1
958    kmx = 1
959    DO k = 1, 360
960      DO i = 2, nlon
961        IF (phy_sst(i,k)<amn) THEN
962          amn = phy_sst(i, k)
963          imn = i
964          kmn = k
965        END IF
966        IF (phy_sst(i,k)>amx) THEN
967          amx = phy_sst(i, k)
968          imx = i
969          kmx = k
970        END IF
971      END DO
972    END DO
973
974    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
975    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
976    ! IM end : verif profil SST: phy_sst
977
978    RETURN
979  END SUBROUTINE profil_sst
980
981END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.